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ABSTRACT 


The search for a universal solution of the equations of motion for a satellite or- 
biting an oblate planet is a subject that has merited great interest because of its 
theoretical implications and practical applications. The discovery of such a solu- 
tion should motivate a reassessment of both the theories that exhibit singularities 
and the physical effects implied by singularities. The practical importance of such 
a solution is the efficiency of simple analytic formulas in predicting simultaneously 
the paths of large numbers of satellites in a multitude of orbits. Here, a complete 
first order solution to the problem of a satellite, perturbed only by the oblateness 
of the Earth, is displayed. The orbit is free of singularities for all parameters and 


is valid for 1000 revolutions with a relative error of the order J? = 10°®. 
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NOTATION 


semi major axis of the initial instantaneous ellipse 
COS to 


eccentricity of the initial instantaneous ellipse (0 < e << 1) 


2cos?6 d¢ 
COB ig at 





constant value of r along an orbit, approximately equal to the 
initial magnitute of angular momentum 

inclination of reference plane 

initial value of 1 


oblateness coefficient of the planet (coefficient of the second harmonic in 


the expansion of the gravitational potential) 

2 (4) 

h?/GM where G is the gravitational constant and M is the mass of the 
planet 

equatorial radius of the planet 

radial distance from the center of the planet to the satellite 

SiN to 

time 

initial value of t 


B 


r 


gravitational potential 

latitude of the satellite 

longitude of the satellite 

angle from line NON’ to the satellite measured in the reference plane 
(Fig. 5.1) 

initial value of 4 

longitude of line NON’ (Fig. 5.1) 

initial value of 


argument of perigee of the initial instantaneous ellipse 
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I. INTRODUCTION 


A characteristic feature of practical satellite orbit prediction is that the engi- 
neer may deal with numerous satellites in a great variety of alternative orbits. Un- 
der these and many other such circumstances analytic relations which can quickly 
approximate an orbit may be far superior to large numerical models. While many 
analytic methods have been developed for the artificial satellite age, most are not 
used in practical orbit prediction because they violate one or more of the following 


principles: 


e The method should provide a solution that is significantly more accurate 
than the two-body solution. 


e The real physical effects of the orbit should be easily distinguishable in the 
solution. 


e The solution should be universal; it should be valid for all orbital parameters. 


The motivation for this research was the desire to develop a method for satellite 
prediction that would embody these characteristics. 

In this analysis, a solution to the equations of motion of a satellite around an 
oblate planet is found by use of a variation of the perturbation technique known 
as the Method of Strained Coordinates. The orbit is valid for 1000 revolutions 
with a relative error of 107°. The solution which is valid for all eccentricities and 
for all inclinations, was obtained by extensive use of the symbolic manipulation 
program MACSYMA. 

The analysis begins with a background discussion of some of the competing 
satellite orbit theories. There is then a development of the equations of motion 


beginning with a derivation of the two-body solution. The various forces which act 


to disturb the two-body orbit are highlighted; a more thorough discussion is given 
for the effects of oblateness. There is a complete treatment of the perturbation 
technique as the equations of motion are solved in detail. The complete first order 
solution is displayed as a function of coordinates and as a function of the orbital 
elements. In addition, an independent analysis of the polar and equatorial orbits 


is performed to serve as a check of the general solution. 


Il. BACKGROUND 


A. INTRODUCTION 


The theory of flight of artificial satellites is closely related to classical celes- 
tial mechanics, one of the oldest and most highly developed branches of science. 
The equations which describe the motion of an artificial satellite are in principle 
identical to the equations of motion of natural celestial bodies. It is not surprising 
then that the results originally derived in classical celestial mechanics have been 
freely used to explain the motion of artificial celestial bodies. 

The foundations of classical celestial mechanics were established in the eigh- 
teenth century when Clairaut, d’Alembert, Laplace, Lagrange, and Euler intro- 
duced theories and analytical methods to explain the large deviation of the Moon 
from an elliptic orbit due to solar attraction. These theories all supplemented 
or complemented the pioneer work that had been done by Newton. Newton had 
correctly indicated the Moon’s variation in eccentricity and inclination and the re- 
gression of the nodes; however, his published theory accounted for only half of the 
motion of perigee. (A century later an unpublished work was found to contain the 
full explanation.) Clairaut in 1749 was on the verge of substituting a new law of 
gravitation for the Newtonian law when he found that second order perturbations 
removed the discrepancy in the motion of perigee. 

In that same century, Euler investigated and developed the perturbative 
function and began the development of the method of variation of parameters, 


which was later extended by Lagrange. 


Lagrange, Laplace, and Poisson all advanced the discipline through their 
investigation of the stability of the solar system. Later, in the nineteenth cen- 
tury, the use of Hamiltonian mechanics was used to great advantage by Delaunay 
whose work influenced the innovators of artificial satellite orbit prediction, namely 
Brouwer, Kozai, and Garfinkel. 

The introduction of the modern computer and the launch of the first artifi- 
cial satellites had a profound effect on the science of celestial mechanics. Classical 
celestial mechanics had been essentially a contemplative science with the principle 
aim being to study the laws of motion of existing heavenly bodies. In contrast, the 
science of the flight of artificial satellites is an active engineering science concerned 
with determining or predicting relatively short term orbits, and in many cases con- 
trolling the satellite’s motion through on-board propulsive devices. It was in the 
exciting climate, following the successful launchings of the first artificial satellites, 


that most of the new methods of satellite orbit prediction were developed. 


B. THE USE OF HAMILTONIAN MECHANICS IN SATELLITE 
ORBIT PREDICTION 


There are certain schools of thought in dynamic astronomy and theoretical 
physics that support the loyal use of Hamiltonian mechanics [Ref. 1, p. 228], 
and many of the new methods in general perturbations take advantage of the 
elegant formalism offered by the Hamiltonian method. The Hamiltonian method 
is referred to here as the formal process of writing the equations of motion for a 


satellite in the canonical form: 














dq, OH dp, —OH 
= = = ] 2 coe S ° 
on i Oe ea) 2) 
where g, = generalized coordinates 


generalized momenta 


3n dq, 
ews eg -E 


1 


Pr 


(H is the Hamiltonian and L is the Lagrangian, L = T — V, kinetic energy - 
potential energy). 
The solution to (2.1) may be written down if a function S can be found, where S 


is any complete solution of the Hamilton-Jacobi equation 


Os OS 
a sre («. st] =0 (2.2) 


It should be noted that (2.2) is tractable only if the variables g, and t are separable 
within S and H. 

Although this present analysis does not use Hamiltonian mechanics to solve 
the equations of motion, a summary of its use is merited since the advances in this 
area have been an invaluable contribution to general perturbation theory. 

In celestial mechanics, one may formulate a Hamiltonian that represents the 
gravitational attraction of the central force and add to it terms which are “per- 
turbing Hamiltonians.” Various sets of canonic variables may be chosen with the 
goal of expressing a zero-order Hamiltonian in a simple form and the higher order 
effects in an iterative fashion. Each term then may be dealt with through a succes- 
sion of canonic transformations. Delauney introduced a systematic procedure for 
isolating parts of the Hamiltonian and then generating a suitable transformation 
in successive steps. A particular feature of his approach is that a periodic term in 
the Hamiltonian may be eliminated with each canonical transformation [Ref. 2]. 
While DeLauney used a procedure that eliminated one term at a time, von Ziepei 
devised a technique that eliminates one angular variable with each transformation. 


This method reduces the number of degrees of freedom while at the same time 


imparting to the transformed Hamiltonian a symmetry of the unperturbed system 
(Ref. 1]. The eliminated variables are referred to as ignorable coordinates since 
they do not participate in the solution of the transformed equations of motion but 
can be recovered after the solution has been obtained. 

Using von Ziepel’s method, Brouwer devised one of the most notable general 
perturbation theories [Ref. 3]. Prior to Brouwer’s method, all previous work in 
artificial satellite theory had written the Hamiltonian as a Fourier series in the 
mean anomaly with coefficients that were infinite series in powers of the eccen- 
tricity. Brouwer used an elliptic approximation for the potential and obtained a 
complete first order theory with some second order development using canonical 
transformations. The essence of Brouwer’s method was to write equations (2.1) as 

dL; OH dl, OH 


a ay dt aL; _ 


where 


[,; =VGMa io ae 
Ly = L1V 1 —- e? i= 


Ww 
Lz = Lecos 1 i..= 0) 
are the Delauney variables where 
a - semimajor axis e - eccentricity 
zt - inclination M - mean anomaly 


w - argument of perigee 1 - longitude of 
ascending node 


Referring back to the Hamilton-Jacobi equation (2.2), the function S is used as 
a generating function to find a new Hamiltonian that leads to simplified canonic 
equations. By choosing S correctly, Brouwer was able to find a canonical transfor- 
mation from the Delauney variables to a set of double primed Delauney variables 


(L;, ’) such that (2.3) have the form 


di dii OH" 
dt dt =«- OL 





(2.4) 


where H* is the Hamiltonian expressed in terms of double primed orbital elements 


by 
L"= J/GMa" I! = M" 
B= Livi-e?® b= w" 
L; = Ly cost" B=" 


The double primed orbital elements are related to the unprimed elements by 


a=a'+ea e=e'+ee 

a tet M=M"+eM 

w=wu'+ew 12=0"4+ E/N 
where the quantities ea,¢e,e1,¢M,ew,eN, are periodic functions of the double 
primed orbital elements. 

Equations (2.4) are solved for the double-primed Delauney variables which 
can then be expressed in terms of the original variables. 

The results initially obtained by Brouwer were not valid at the critical incli- 
nation of 63°.4, and they were questionable when either the inclination or eccen- 
tricity were near zero [Ref. 4]. O. K. Smith devised a method for dealing with the 
problem of zero inclination and eccentricity [Ref. 5], but Brouwer subsequently 
challenged the validity of his method. Later, Lyddane |Ref. 6] was able to suc- 
cessfully remove the restrictions on small values of eccentricity and inclination by 
reformulating Brouwer’s work in terms of an alternate set of variables. 

Brouwer used the central force term as the first approximation for the po- 


tential since there is no exact solution to the equations of motion of a satellite 


under the influence of the more complex potential described by an oblate planet 
having axial symmetry. Other authors have attempted to introduce a new poten- 
tial which approximates the Earth’s potential better than the central force term 
alone and also leads to an exact solution. Most notable has been the work of 
Sterne (Ref. 7] and Garfinkel [Ref. 8]. Sterne’s potential function accounts for 
most of the standard potential through the second harmonic and it leads to a 
solution of canonical constants that are free of first order secular perturbations. 
The remaining effects of the earth’s oblateness and other forces are allowed for 
in the perturbing Hamiltonian which causes the six canonical constants of the 
unperturbed solution to undergo variations with time. Sterne provided the inspi- 
ration for Garfinkel’s method which is essentially the same as Sterne’s but more 
developed. Garfinkel included the second and fourth harmonics and arrived at 
a solution that is reducible to quadratures. Garfinkel’s original solution was not 
valid at the critical inclination, but in a later paper he removed the singularity 
through a variation on his perturbation theory |Ref. 9]. 

While Garfinkel did his analysis in spherical coordinates, Vinti (Ref. 10], 
derived a potential expressible in oblate spheroidal coordinates. In his original 
analysis, Vinti introduced a potential function and associated coordinate system 
that would lead to the separability of the Hamilton-Jacobi equations. In a subse- 
quent work |Ref. 11], Vinti showed that the equations result in a solution reducible 
to quadratures. The Vinti potential is an exact expression of the Earth’s potential 
through the second harmonic. The theory provides perturbed coordinates, not 
perturbed elements. In his second order analysis [Ref. 12], Kozai criticized this 
characteristic, and he chose not to use the Vinti potential since it would require 


changing the definitions of the conventional orbital elements. 


Morrison |Ref. 13] showed that the von Zeipel method is a particular case of 
the method of averaging, and Liu [Ref. 14] used the latter technique to study the 
combined effects of air drag and planet oblateness on a satellite orbit. The method 
of averaging, unlike von Ziepel’s method, does not require that transformations 
be canonical. The method of averaging has been used extensively in recent years 
by Lorell and Liu [Ref. 15], McClain [Ref. 16], and Hoots [Ref. 17]. However, the 
validity of the method has been challenged; most notably Taff [Ref. 18] doubts 
its rigorous foundation. Arnold [Ref. 19] notes that the principle of averaging is a 
vaguely formulated and rigorously untrue assertion, but he adds that sometimes 
such assertions are fruitful mathematical sources. It should be also noted that 
the solutions obtained in Liu’s analysis [Ref. 14] are not valid for near circular 
equatorial orbits nor at the critical inclination, while those obtained by Hoots are 
valid only for small eccentricities (0 < e < .1), and they are invalid at inclinations 
near 0° or the critical. 

Hamiltonian mechanics has provided a rich source of literature in orbit the- 
ory; however its practical applicability has been questioned in some textbooks. 
Roy |Ref. 20] briefly outlines the use of the canonic equations, while Taff [Ref. 
18p. 322] states that he does not see any additional practical applications provided. 
Baker |Ref. 21] chooses not to represent the subject. A general criticism is that 
the process of generating suitable transformations in the perturbation procedure 
tends to make the coordinates and the momenta less distinguishable on physical 
grounds and more difficult to relate to the set of natural coordinates which were 
used to write down the initial set of differential equations. While the ultimate 
form of the governing equations may be simple to solve, there remains the tedious 
task of obtaining explicit results in terms of physically meaningful coordinates or 


elements. 


C. KOZAIY’S METHOD 

Prior to his work cited above, Kozai (Ref. 22] developed a method for find- 
ing the perturbations on the orbital elements of a satellite considering only the 
oblateness of the Earth. Kozai developed a disturbing function based on the 
Earth’s departure from a sphere, and he used a version of Lagrange’s planetary 
equations to formulate the solution. Kozai used the standard form for the Earth’s 
potential and included the harmonics J2, Js, and J,. Despite the use of the higher 
harmonics, the theory is first order. Kozai expressed short-period terms in J2, the 
secular in Jz, J,, and J?, and the long period terms in J2, ie and ae The analytic 
expressions are developed using the standard orbital elements. 

Kozai’s work is cited here because, due to its simplicity, the method has be- 
come very popular in many textbooks and handbooks on orbital mechanics. [Refs. 
20, 23, 24]. However, Taff cites Kozai’s method as an example of misapplication 
of perturbation theory [Ref. 18p. 332]. Taff challenges the assumptions made by 
Kozai in his analysis, and he points out that the method is invalid at the critical 
inclination. 

As was stated in Chapter J, a motivation for this current analysis was the 
purpose of finding a perturbation method that would lead to universal solution. All 
of the methods discussed thus far have particular problems at certain inclinations 
or eccentricities. Some of the problem cases have been resolved by unique efforts 
(These cases are discussed in Appendix C.); however, one should question the 
underlying validity of any perturbation method that produces singularities. There 
has been no satisfactory way found for avoiding the critical inclination singularity 


in Kozai’s method. 
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D. THE DIRECT METHODS 

R. E. Roberson [Ref. 25] devised an approach for finding the qualitative 
and approximate quantitative results concerning the behavior of a set of orbital 
elements in the gravitational field of an oblate planet. The motion of a near 
satellite around the planet is simply described to first order by introducing a frame 
of reference which contains a mean orbital plane having a constant inclination. 
Both the reference frame and the orbit plane rotate at a constant angular velocity. 
With respect to this doubly moving reference frame, the motion of the satellite 
differs from pure elliptic motions by only periodic perturbations. 

King-Hele [Ref. 26] advanced the approach taken by Roberson. He intro- 
duced a non-rotating reference frame with a orbit plane that continually rotates 
at a non-constant rate about the Earth’s axis. A relation is found between the 
rotating orbit plane and the angular rate of travel of the satellite. The equations 
of motion are written in position coordinates, but are subsequently rearranged 
in terms of a modified set of orbital elements. The inclination of the rotating 
reference plane is held strictly constant in the analysis. King-Hele formulated 
the problem by employing a power series expansion in terms of the eccentricity; 
therefore, the method is limited to small eccentricities. The final solution con- 
tains an incomplete set of workable elements and hence the method gives mostly 
a qualitative description of the satellite’s behavior due to oblateness. 

King-Hele’s analysis was the inspiration for the work of Brenner and Latta 
(Ref. 27]. They improved King-Hele’s original analysis by abandoning the condi- 
tion of constant orbit inclination and by retaining the eccentricity in closed form 
expressions. Brenner and Latta limited their analysis to small eccentricity al- 
though the method is not so restricted. They obtained an approximate first order 


solution and demonstrated that the method is valid for higher order analysis. 
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An advantage of the direct method is that one may use ordinary perturbation 
analysis (the Method of Strained Coordinates) to solve the equations of motion. In 
addition the chosen set of orbital elements used throughout the analysis correspond 
closely with the classical elements; therefore, physical interpretation is facilitated. 

The disadvantage in using the method is shared by most all other procedures 
that must use a perturbation scheme. The process requires the manipulation of 
massive algebraic expressions. However, this drawback has been greatly reduced 
by the introduction recently of large symbolic mathematics programs, such as 


MACSYMA, that handle the bookkeeping. 
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lil. THE TWO-BODY PROBLEM 


A. INTRODUCTION 


The central problem of celestial mechanics is the two-body (or Kepler) prob- 
lem. Simply stated, the problem is to solve for the motion of two particles inter- 
acting through their mutual gravitation in an isolated space. A solution of the two 
body problem often represents physical reality in an acceptable way. For instance, 
the orbit of the Earth around the Sun may be treated, to a first approximation, as 
a two-body problem because the influence of perturbing bodies, the Moon, Jupiter, 
etc., are small compared with the Sun’s gravitational attraction. Likewise, more 
complex problems such as a spacecraft mission to Mars, a four body problem — 
Earth, Sun, Mars, spacecraft, may be treated by breaking the flight into three 
two-body problems. Such a technique is used in the patched conic method where 
two-body solutions are literally patched together. 

There are other more important reasons for studying the two-body problem. 
It is the only gravitational problem in dynamics, other than very specialized cases 
in the three-body problem, for which there is a complete and general solution, 
and it is possible to gain considerable insight into more general phenomena of 
motion by a thorough study of the two-body problem. In fact, the most complete 
theories of celestial motion use functions appearing in the solution of the two- 
body (elliptic case) problem as elementary functions. The solution is central to 
this present analysis because it will serve as a starting point for the generating of 
analytical solutions that are valid to higher orders of accuracy. These solutions, 


called general perturbation theories, are the subject of Chaper VI. 
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In this chapter, the equations for the two-body problem will be derived and 
solved. The first step is to chose a coordinate system in which the laws of Newton 
hold (an inertial coordinate system). In practice, the reference frame of the “fixed” 
stars provides a very good approximation to an inertial reference frame. Next, 
following the method of Nelson and Loft [Ref.28 pp.82-84] it will be shown that 
the center of mass of two bodies in this coordinate system is also an inertial 
reference point. Then the equations of motion will be derived and solved using 


the polar angle @ as the independent variable. 


B. THE DIFFERENTIAL EQUATION 

Figure 3.1 shows two mass centers at position r; and rz. The origin O is 
defined to be an inertial reference point. The distance between the two mass 
centers is r where 


r= fre —TPy. 


Combining Newton’s third law of motion and his law of universal gravitation gives 





Figure 3.1: Two Body Problem 
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the force equations for the two bodies: 


Tye S Gm, mor 
dt? rs 
d’*r. —Gm mz r 
at - a 


where G is the universal gravitational constant. 
Adding (3.1) and (3.2) and integrating results in 
dr, dr» 


m,— + m,— = constant. (3.3) 


dt dt 


Now the center of mass of the system (barycenter) in Fig. 3.1 is defined as: 


(my ae m2)¥o = mr, + MPL. (3.4) 
By combining equations (3.3) and (3.4) the following result is obtained, 
dro my dry m2 2 
— = ———_ —— + ————_.-——~ = constant. 
dt m, +m, dt m, + my, dt 


So the barycenter moves with constant velocity, and it too is an inertial reference 
point. Subtracting (3.1) from (3.2) results in 
dr, d’r, dr —G(m,+m,)r 


di? ae gi 3 (9 








the solution of which gives the position of either body relative to the other. Now 


choose the barycenter as the origin and define the position of m, and mz, respec- 


tively, as 
rol = <1 a Ome 
mm, + mM 
my 
loo = 1) — Toe 
my, + M2 
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Substitution of the above equations in (3.1) and (3.2) yield separate equations for 


the motion of each body relative to the barycenter: 











dro d’roy _ d*ro, = —Gm} ~ 
dt? ' dt? dt?  (m,+m,)?r3, ° 
(3.6) 
dro, d’roz_ @'ton _ Gm 
dt? ' dt? dt? ~~ (m,+m,)?r38, 


Equation (3.5) and equations (3.6) are of identical form, differing only by a con- 


stant, so that 


2 
a + — =0 (3.7) 
is the vector differential equation of motion for either of the two bodies. r is the 
distance to the other body, or to the barycenter, according to the appropriate 
choice of GM. For the problem of a satellite orbiting the earth, the mass of the 


satellite can be neglected in comparison with the mass of the earth, therefore GM 


is the product of the universal constant and the mass of the earth. 


C. THE INTEGRATION OF THE TWO-BODY PROBLEM 
There are no cross products involved in equation (3.7); therefore, all motion 
must lie in the plane that contains r and x The scalar components of acceleration 


are 





d*6 m 2dr dé : an 
l ee oe 
dt? dt dt aa 
dr d)\*  —GM a 
dt? \Gi) ag) 
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(r? a) = O and integrating yields 


Writing (3.8) as ig dt 


dé 
r?-— — h = constant (3.10) 


dt 


where A is the specific angular momentum. 


Equation (3.10) is an exact integral of (3.8). It corresponds to Kepler’s 


empirical law of constant areal velocity which states that the area swept out by 


the radius vector of a planet is uniform in time. 


From (3.10), the independent variable may be changed from t to 8, e.g., 








a@_ha 
dt r?2d6 
o (3.9) becomes 
ld’r | 2 (ar\" 1 _ 6M _@ 1.) ,1_ GM 
r?2 de?" +3 \ dé ~ Re gee TT 


Since this equation is linear in the reciprocal of r, it may be written as 


d*u GM 
ap? are 


which has the solution 


GM 
errant Acos(@ — w) ell) 





or reintroducing r, (3.11) becomes 


h?/GM 


~ 14+ (Ah?/GM) cos(6 — w) 


which may be written as the polar equation of a conic section 


Las, 


p 


rr cos(@ — w) 28 


so that 


p = h?/GM and e = Ah’/GM. 


D. ELLIPTIC MOTION 

Equation (3.12) is the equation of a conic with prime focus at O. The conic 
has a semi-latus rectum h?/GM, an eccentricity e, and a semi-major axis a that 
makes an angle f = 6 — w with the horizontal axis (Figure 3.2). The extreme 
endpoints of the major axis of the orbit are referred to as apsides or apses. The 
point nearest the prime focus is called perigee and is given by 6 = w. The point 
farthest from the prime focus is given by 6 = w + 180° and is called apogee. The 
angle w, “argument of perigee”, will be discussed later in this chapter. 


The energy of the satellite in the orbit is conserved and is equal to 


E=m,(v?/2—GM/r) =m,C 








- 7a 
Origin at focus: 






Figure 3.2: The Elliptic Orbit 
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where C is the total energy of the satellite, v?/2 the kinetic energy and —-GM/r 
the potential energy of the satellite, all per unit mass (m,). In general C is equal 
to —-GM/2a, so that the satellite’s energy depends only on the semi-major axis. 


From this relationship it is easily shown that the satellite’s velocity is given by 


v’ = GM(2/r — 1/a). 


Then the relationships r, = a(1 — e) and r, = a(1 +e) result in 


co GM(1 +e) 
‘ a(1 —e) 
and 
nL. GM (1 — e) 
a(1 +e) 


so that the velocity is a maximum at perigee and a minimum at apogee. 


The area of the ellipse is ma’,/(1 — e?), and the rate of description of area is 


rae = 7 Since h? = GMa(1 — e’), the orbital period may be written as 





By defining the mean motion n as T = %, so that n?a* = GM, one may proceed 
to derive an expression for position versus time in the elliptic orbit. 

The orbital ellipse APB with center at C touches at perigee, A, and at 
apogee, B, a concentric circle also centered at C' which has as a radius the semi- 


major axis of the ellipse. The circle C is known as the auxiliary circle, and is 


geometrically related to the ellipse by the relation 


PN = P'N,/(1 — e?) 
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where e is the eccentricity, P any point on the ellipse, N the foot of the perpendic- 

ular through P upon AB, and FP’ the intersection of this perpendicular with the 

circle C. The angle ACP’ = E is known as the eccentric anomaly (Figure 3.3). 
Let 7 be the time of perigee passage and ¢ the time, then t —7 is the time 


since perigee passage. The quantity, 


M =n(t—7), (3.13) 


is called the mean anomaly. Using the geometry of Figure 3.3, one may derive the 


following relationship between the anomalies: 
M = E -esin E, (3.14) 


which is known as Kepler’s equation. If the position of a satellite in a fixed orbit 
relative to the earth is desired at some specified time t, then equation (3.13) gives 


M and (3.14) gives E . The distance to the satellite is found by the relationship 


r=a(l-—ecos E). 


The angular position of the satellite is defined by the true anomaly, 6 —w = f,or 


the angular distance from perigee: 





Figure 3.3: The Eccentric Anomaly 
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fan(j7/2)— eae ta & 
In actual practice, the second step in the above process, solving (3.14) is a bit 
more involved since a closed form solution to Kepler’s equation does not exist. 
However, dozens of methods of successive approximations have been devised. For 


instance, by use of a numerical method M can be calculated for a few values of E 


and then the correct E can be found by inverse interpolation. 


E. CONSTANTS OF THE MOTION 

In section B, the original equations of motion (3.1) and (3.2) were reduced 
o (3.7). Thus the problem was reduced from one of three second order differential 
equations requiring twelve constants to one of three second order equations with six 
constants. A discussion of the constants of motion, some of which were introduced 
during the solution of the two-body equation, is the purpose of this section. The 
six constants may be written in a variety of forms, a choice among forms is usually 
made with the purpose of simplifying the problem. 

Equation (3.7) was solved by the classic technique of changing the indepen- 
dent variable from t to 6. But (3.7) in its present form is integrable; that is, 
there exists sufficient time independent first integrals, or functions that are con- 
stant along the motion, to specify each orbit. The first of these is the angular 


momentum. Cross-multiplying (3.7) by r results in 


Z 


Ee a 
or 
ez 0 
r —_ = 
dt? 


za 


And since 


d dr dr dr d*r 
& (rx $ See og, RANT 
then 
<r xv) = 0. 
By integration 
YX — 


where h, the angular momentum, provides three constants of motion. Similarly, 


by cross-multiplying (3.7) by h, one may obtain the Lenz vector 


ef peta ($= =o) 
dt GM fr dt 

where e is a vector along the major axis of the orbit pointing toward the position of 
closest approach or perigee (|e] = e). The vector e provides only two independent 
constants since h and e are perpendicular vectors (h-e = 0), and one remaining 
constant is required. 

The vector integrals h and e specify that the orbit will lie in the plane per- 
pendicular to h and have a shape determined by e. The classical orbital elements: 
a (semi-major axis), e (eccentricity), ¢ (inclination), N (longitude of the ascending 


node), and w (argument of perigee), may be derived from these vectors. 


From (3.12) the equation for r may be written as 


h?/GM 
p= 
1+ecos f 


where f is the angle between e and r (the true anomaly). 


Restricting e to the elliptic case results in 
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h? h? 
= ir sv ——— 
GM(1+e)~ ~ GM(1-e) 


so that the major axis of the ellipse is 


2h? 
GM(1 — e?) 


20.— 5) ro = 


Knowing e and a gives the shape and size of the orbit. The orbit now can be 
oriented in a coordinate system. Reference is made to Figure 3.4 and the two 
angles : and (2. : is the inclination of the orbit plane defined as the angle between 
the equatorial plane and the orbit plane. Since this is the same angle as the angle 
between the z axis (k unit vector) and the angular momentum vector h, 1 may be 


found by 


cost= —. 





Figure 3.4: Constants of the Orbit 


23 


The angle , the longitude of the ascending node, is the angle in the equato- 
rial plane, between the i unit vector and the point N where the satellite crosses 
through the equatorial plane in a northerly direction measured counterclockwise 
when viewed from the north side of the equatorial plane. The vector n lies along 
NO such that 

n=k xh. 


Therefore M may be found by 
tl. 
cos = —. 
n 
As was stated above, the vector e, points toward perigee. The angle w (the 
argument of perigee) measures the distance between NO and perigee and can be 
found by 
n-e 


cos W = —. 
me 


With the constants a,e,1,2,w specified, the orbit is defined in the coordinate 
system. The remaining task, to locate the satellite in the orbit at any time, 
requires one more constant. 

The final constant of motion is given by the relationship between the mag- 


nitude of the angular momentum and the true anomaly (f): 


t , 
[ na= rdf. 


By a change of variables from f to the eccentric anomaly EF, the above equation 
may be easily integrated. The result by this analytical method is identical to 
equation (3.14), Kepler’s equation, which was previously derived by geometry. 


The constant 7, time of perigee passage, is the final constant of integration. 
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F. SUMMARY 

This completes the analysis of the two-body problem. It has been shown that 
a combination of six constants will strictly define the motion of a satellite under 
the influence of a central gravitational force. The six that were chosen, referred 
to now as the orbital elements, can be used to find other constants of the motion 


including r and v, the distance and velocity vectors of the satellite. 
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IV. SATELLITE PERTURBATIONS 


A. INTRODUCTION 

As was demonstrated in the last chapter, the classical two-body problem has 
solutions that can be written in closed form when the polar angle (or the eccen- 
tric or true anomaly) is used as the independent variable. If an additional force 
acting on either of the two bodies is introduced, the resulting equations of motion 
usually no longer have closed-form solutions. When the magnitude of such a force 
is small compared to the central gravity term, the force is called a perturbation. 
The resulting orbit experiences a small departure from the Keplerian orbit, at 
least initially. These departures are also called perturbations. Under certain cir- 
cumstances, it is possible to make analytic approximations to the effects of the 
perturbing forces, though a precise solution cannot be obtained. Generally, the 
methods consist in determining the exact equations of motion and then assuming 
that their solutions do not depart appreciably from the case of no disturbing force. 
Then only an indication of the actual motion of the body can be obtained. Pre- 
cise solutions can be found for specific initial conditions by numerical integration 
techniques, but these solutions give little insight into the dependence of the mo- 
tion on the parameters of the disturbing force. In some cases, the approximations 
obtained with analytic methods may exceed the precision of numerical methods, 
especially if the prediction is required for a long period of time and there is a clear 


dominance of one particular force. 
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In the case of a close satellite about a non-spherical planet, a potential func- 


tion V can be formed such that 
VS Vo te R 


where Vo is the potential function due to the two-body problem and R the dis- 
turbing function that is at least an order of magnitude less than Vo. Many general 
perturbation theories make use of the fact that the two-body orbit of the body due 
to Vo changes slowly due to R, and they attempt to obtain analytical expressions 
for the changes in the orbital elements due to R within a specified time interval. 
If the elements of an elliptical orbit are ao, €o,19,Mo, Wo, and to, the ellipse with 
these elements is referred to as the osculating elements at tp. The velocity of the 
disturbed satellite at this time in its osculating ellipse is equal to its velocity in 
the actual orbit. 

At a future time t,, the elements will change due to the presence of R. 
For instance, a) will be changed to a;, and the quantities (a; — ao), (e€1 — €0), 
etc. are called perturbations in the elements. Corresponding to the perturbations 
in the elements are perturbations in the coordinates and velocity components. 
There are, however, at least two reasons for using orbital elements rather than 
coordinates to describe the motion of the satellite. First, the elements do not 
exhibit a variability of anomalistic motion that the coordinates do, therefore any 
variation can be attributed directly to the perturbing forces. Second, the elements 
give a clearer geometric picture than do the coordinates; hence the effect of the 
perturbation on the orbit can be seen immediately. 

There are various kinds of disturbances that an orbit can experience, the 


severity of each is usually due to the altitude of the satellite. It is the purpose of 


a 


this chapter to give a qualitative description of the most important disturbances, 


and to relate the relative magnitudes of each. 


B. THE EARTH’S GRAVITATIONAL FIELD 

The two-body problem assumes that the earth is a sphere, while in reality the 
earth is flattened somewhat at the poles and bulges correspondingly at the equator. 
Such a shape is called an oblate spheroid. In the science of geodesy, it has been 
useful to define a reference ellipsoid as a mathematical surface which is an idealized 
approximation to the earth’s actual surface. The study of satellite orbits has 
established a flattening of the terrestrial ellipsoid as 1/298.24, which corresponds 
to a difference between the equatorial and polar radii of 21.4 kilometers. 

Another surface commonly used in geodesy is the geoid, which is the equipo- 
tential surface that coincides on the average with mean sea level in the oceans. The 
geoid is everywhere perpendicular to a plumb-line since gravity is always normal to 
its surface. Before the advent of artificial satellites, it was generally accepted that 
except for relatively small variations resulting from the presence of mountains or 
deep valleys, the geoid could be regarded as approximately an ellipsoid. Now, the 
surface of constant gravitation can be more accurately portrayed by representing 
the potential as a series of quantities known as spherical harmonics, each of which 
makes a contribution, positive, negative, or zero, to the total. The contribution 
of any harmonic is determined by a factor, usually represented by the symbol J 
and commonly referred to as the value of that harmonic. These J values for a 
planet’s gravitational field can be determined from observations of a satellite orbit 
and they can be related to the shape of the geoid. 

A large number of harmonics may be required to precisely represent a planet’s 


gravitational field, but in practice the higher harmonics make such a small contribution 
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that they can be neglected, at least to a first approximation. The zonal harmonic 
Jo expresses the overall size of the geoid, while J,, the first degree harmonic de- 
termines the center point of the geoid in the north-south direction. The other 
harmonics represent deviations from the spherical shape as shown by Figure 4.1. 
It is seen that the contributions from the even harmonics are symmetrical about 
the equator, while the odd harmonics corresponds to contributions that are asym- 
metrical. The degree of the harmonic gives the number of undulations in the shape 


of the surface. 


- 





~~ 


Figure 4.1; Qualitative Representation Of The Harmonics Of The Geoid 


The results given thus far consider only the zonal harmonics, which are in- 
dependent of longitude. The tesseral harmonics give east-west deviations from 
symmetry. Satellite observations of the tesseral harmonics have led to the con- 
clusion that the equator of the earth’s geoid is slightly elliptical, the difference 
between the longest and shortest axes being about 400 meters. 

The following tables from Kozai [Ref. 29] gives representative values for the 


coefficients of the earth’s zonal and tesseral harmonics. 
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TABLE 4.1: ZONAL HARMONICS 


n J, X 10° n J, X 10° 

2 1082.48 + . i 
3 —2.57 + 
4 
5 


—1.84 + 
—-.06 + 





TABLE 4.2: TESSERAL HARMONICS 
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s 
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m hm Bm CO CO WF NO 
WON BS WN ee 





C. ATMOSPHERIC DRAG 


Atmospheric drag is the most complex and the most difficult of the artificial 
satellite perturbations. The complexity arises from the fact that an exact force 
law is not known and the atmosphere is variable. This variability results from 
the fact that the atmosphere is not spherically symmetric and that lunisolar tides, 
diurnal heating, strength of the solar wind, etc., all effect atmospheric density. 


In addition, the atmosphere is moving with the rotating earth. Therefore this 


perturbation is difficult to model analytically or numerically. 


For most Earth satellites, atmospheric drag removes the satellite’s energy, 
thereby causing it to drop in altitude and thus increase its speed. The increased 


speed and the lower altitude increases the drag force with the result that the 


satellite spirals into Earth. 
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Satellite drag is measured along the same direction as the velocity vector of 


the satellite. The formula used is 


= 
Fp = > CoAely —Val(Va —V) 


where v is the satellite’s inertial velocity, vg is the inertial velocity of the atmo- 
sphere, p is the atmospheric density, A is the effective cross-sectional area of the 


satellite, and Cp is the drag coefficient, Cp ~ 2. 


D. SOLAR RADIATION PRESSURE 

In contrast to drag which is a tangential force, solar radiation pressure is 
radial. By solar radiation on a satellite is meant the net acceleration resulting 
from the interaction (i.e., absorption, reflection, or diffusion) of the sunlight with 
the surface of the satellite. In general, the illumination of a low altitude satellite 
orbit will be non-uniform due to the Earth’s shadow, albedo, and atmospheric 


absorption. The magnitude of the solar radiation pressure is given by 


_AP 
4mmcD? 


where A is the effective cross-sectional area of the satellite, P is the total radiated 
solar power, m is the satellite’s mass, c is the speed of light, and D is the satellite- 
Sun distance. 

It can be seen that the perturbing effect of solar radiation particularly effects 
satellites with a large area to mass ratio. Such was the case of the satellite Echo I, 
which when inflated was a sphere with a total exposed area of about 31,000 square 
feet and weighed only 135 pounds. Its initial orbit was approximately circular with 
an altitude of 1000 miles, but the pressure of solar radiation brought down the 


perigee to 600 miles at times. 
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A. OTHER PERTURBING FORCES 

Gravitational perturbations due to other celestial bodies (the Moon, the Sun, 
and the planets) are caused by the differences between the force on the Earth and 
that on the satellite or the “tidal” force. Perturbing forces from the Sun and Moon 
become significant as the altitude of the satellite increases. Planetary perturba- 
tions are very small, with Venus and Jupiter providing the largest contributions. 

Since classical mechanics is used in this analysis, relativistic corrections may 
be regarded as a small perturbation to the motion. Other minor perturbing forces 
include atmospheric lift and electromagnetic forces. 

The following table from Milani [Ref 30] gives the magnitudes of various 
perturbations on three different satellites: 


ACCELERATIONS ON SPACECRAFT AT VARIOUS 
ALTITUDES (cm/sec?) 


TABLE 4.3: 


LAGEOS SEASAT 





Geosynchronous 







satellite 
a= 42,160KM 12,270 7,100 
A/M = lcm g .007 2 









Earth’s >a iG, 7.9 x 10° 


monopole 


2.8 x 10 





















albedo 





Earth’s t4SoOloe: 1.0% 10> 9.3 alos: 
oblateness 
Higher 4 Oat 8.8x10°© 7.0 x 1074 
harmonics 
Moon 1S Se 1 Fie Lo oclOes 
Sun 3.3 x 1074 96x107> 56x 1075 
Venus 4.3 x 107 1S lee oe) Sa vaaliOu” 
General Zo Oe 93x107% 4.9x 107’ 
Relativity 
Air Drag 0 (7) 3 Oparee 2.0 hee 
Solar 4.6 x 107° 3.2 103 9.2 10 
Radiation 
Earth’s A 10.” 3.4x107% 3.0x 107° 





V. DEVELOPMENT OF THE EQUATIONS OF 
MOTION 


A. PRELIMINARIES 
1. Overview 
In this chapter, the groundwork will be laid for solving the equations of 
motion for a satellite about an oblate planet. To begin, a discussion of the assump- 
tions made in the analysis will be given. Second, the special coordinate system, 
which was introduced by King-Hele [Ref. 26] and refined by Brenner and Latta 
[Ref. 27] will be developed, followed by a derivation of the relationships among 
the astronomical angles of the coordinate system and the angles of a spherical 
coordinate system. Next, an expansion for the potential of a planet modeled as 
an oblate spheroid will be derived. Finally, the equations will be transformed so 
that they can be solved by an ordinary perturbation method. 
2. Assumptions and Limitations 
The basis for the assumptions made in this analysis is the requirement 
that a solution to the equations of motion for a satellite orbiting a planet be 
accurate to within a relative error < 10~°. Therefore, all perturbative forces 
that are of magnitude 10~° or smaller compared to 1 may be neglected. This 
requirement then allows one to model the earth as an oblate spheroid with an 
axially symmetric gravitational potential. This is to say that all zonal harmonics 
except the second in the expanded potential may be neglected. In addition, all 


coefficients of the tesseral harmonics are small enough to neglect. 
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The assumption of an axially symmetric potential is a reasonable as- 
sumption if the earth is used as an example of an oblate planet. As was suggested 
in Chapter IV, if one normalizes the expanded geopotential of the earth such that 
the contribution from the central force (1/r?) is 1, then the contribution from the 
zonal harmonic for the oblateness contribution is ~ 107°, and all other harmonics 
(zonal and tesseral) are < 10°. 

The gravitational effects of the Sun and Moon are also neglected. This 
assumption will remain valid out to a distance of at least 6,000 kilometers above 
the Earth, where the oblateness perturbation is still about 1000 times larger than 
the attraction of the Sun. However, at altitudes near geosynchronous (35,800 km) 
the perturbative forces of the Sun and Moon are nearly equal to that of oblateness, 
and therefore would have to be considered. 

The most important limitation to the analysis is the neglect of atmo- 
spheric drag. For high altitude satellites with small eccentricity, the neglect of air 
drag is unimportant; however, for very low altitude satellites the effect of air drag 
would begin to dominate the oblateness force after a number of revolutions. Also, 
high altitude satellites with highly eccentric orbits would be greatly affected by 
drag as they pass through perigee. For these particular cases, the desired accuracy 
for 1000 revolutions would not be achievable; however, in many cases the solution 
could be accurate for a few orbits. As was discussed in Chapter IV, air drag re- 
mains the most difficult perturbation to model. Since an accurate geopotential 
model allows for a better determination of fluctuations in drag of the atmosphere 
through which the satellite is moving, the analysis conducted here is valuable even 
when air drag becomes the dominating factor. 

Also neglected in this analysis, are the remaining perturbative effects 


mentioned in Chapter IV. Solar radiation pressure, the Earth’s albedo radiation 
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pressure, and the general relativistic correction to Newton’s law of gravity are all 
neglected since their relative contributions are small. 

The final limitation to the analysis is the assumption that the differential 
equations which describe the motion of the satellite have a solution asymptotic to 
the form specified in the analysis. As long as the true solution does not depart 
from the form specified, the analysis is valid. 

3. Order of Magnitudes 

Throughout this analysis, reference will be made to the relative mag- 
nitudes of the perturbing forces; hence it is important to establish a convention 
for orders of magnitudes. The approach established by King-Hele [Ref 26] will be 
followed here with the exception that a small eccentricity is not assumed. Here, 
J = Jo, the coefficient of the second zonal harmonic in the potential is used as the 


basic small unit. Mathematically, 


f(J,8) = O(776") 


means that there exists for all 


Or <= 1Omeandio = 6 (2710: 


a K independent of J and @ such that 


| Fleet Or 
The following terms are defined: 
Zero order = O(1) ~ 1.0 
First order = O(J),O(J?6),O(J°67),... ~ 1073 
Second order = O(J*),O(J°0),O(J‘67),... ~ 107° 
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B. THE COORDINATE SYSTEM 

Figure 5.1 is an inertial reference system of spherical coordinates with the 
prime direction pointing toward the vernal equinox at epoch 1900.0. The equato- 
rial plane, latitude (5), and longitude (¢), as shown in the figure have their usual 


meanings. 


reference plane 





equatorial plane 





Figure 5.1: The Reference System 


As was demonstrated in Chapter III, the path of a satellite of a strictly spherical 
planet lies entirely in a fixed plane, and the motion of the satellite is described by 
the solution to the two-body problem. With the angular momentum vector fixed 
in space, the intersection of the orbit plane and the equatorial plane describe the 
fixed angle measured from the prime direction to the intersection (or node). 
The longitude of the ascending node (NM), the inclination (1), and the argument of 


perigee (uw), fix the orbit plane in the coordinate system. 
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The effect of the oblateness perturbation is, in general, to move the satellite 
out of the original two-body orbital plane. One specific effect is to rotate the 
orbit plane about the planet in the opposite direction to the satellite’s motion. 
The physics of this motion is easily demonstrated in Figure 5.2. Oblateness is 
represented by an additional mass about the equator of the planet. The radius 
vector r and the satellite S lie in a plane which is named here the “reference 
plane”. If there were no equatorial bulge, the direction of the gravitational force 
would coincide with the radius vector and the angular momentum vector would 
remain in a constant direction normal to the plane. The equatorial bulge; however, 
adds a component of force that does not lie along r. This additional force adds a 
torque r =r x F. The direction of r is into the paper at S for a prograde orbit. 
Therefore as the satellite orbits, the angular momentum vector rotates about the 
Z axis. The reference plane also rotates, and the rotation is measured by the 


change in 9%. 






Equotoria! 
plone 


\ Oblate 
spheroid 


Figure 5.2: Rotation of the Reference Plane 


The line NON' in Figure 5.1 is referred to as the “line of nodes”, and the 


ascending node describes the point where the satellite enters the northern hemi- 


3/ 


sphere. For the two-body problem with a non-rotating plane, the line of nodes is 
a fixed reference line. Now, for the rotating plane, when the satellite is actually 
at a node, the line of nodes passes through the satellite, but at other times the 
definition “line of nodes” is an arbitrary one, and the angle 2 is simply given as 
a function of time. 

Another motion of the reference plane that is a result of the oblateness 
perturbation is an in-plane rotation called precession. As was demonstrated in 
Chapter III, a particle executing bounded, non-circular motion in a central force 
field will always have a radial distance from the force center that is bounded by 
Tmazr <1 <Pfmin» Lhat is, r is bounded by the apsidal distance. Such an orbit is 
called a closed orbit and it is characterized by the fact that all apsidal angles are 
equal. For instance, the apsidal angle is 7 for elliptic motion. But if the radial 
dependence deviates slightly from 1/r’, then the apsides will precess or rotate 
slowly in the plane of motion (Figure 5.3). This motion is analogous to the slow 
rotation of elliptic motion of a two dimensional harmonic oscillator whose natural 
frequencies for each dimension are almost equal. The rotation of the line of apsides 
is also known as the precession of perigee since the value for @ for which r is a 


minimum varies in time as the apsides rotate. 


Wwe eo 


~~ 





Figure 5.3: Rotation of the Line of Apsides 
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For a satellite orbiting an oblate planet the rotation of the line of apsides is 
opposite the direction of satellite motion for inclinations less than sin?! /4/5 and 
in the same direction as the satellite for inclinations greater than sin™! (4/5. If 
the inclination is exactly sin™! /4/5 then there is no rotation of the line of apsides, 
and perigee remains constant. This inclination is known in Astrodynamics as the 
“critical inclination”. A detailed discussion of the critical inclination is contained 
in Appendix C. 

A final motion of the reference plane caused by the oblateness perturbation 


is the periodic variation of the inclination about the initial inclination 1p. 


C. ANGLE RELATIONS 

As a preliminary to deriving the equations of motion, it is necessary that 
the relationships among the spherical coordinates and the angles describing the 
reference plane be established. 

Reference is made to Figure 5.4, where i,i’,j, and j’ define the equatorial 
plane. Let a,b, and c denote three unit vectors, where b and c are in the reference 
plane. c points to the initial point @ = 2/2 in the reference plane where @ is 
measured from the line of nodes. b points to the position of the satellite, S. a is 
perpendicular to the reference plane and therefore points in the same direction as 
the angular momentum vector h. a and c are both perpendicular to the line of 


nodes. Then 


b=cos ¢cosé6i+sing cosé6j+sinék. 


Now measure b from the line of nodes. 


b = cos(¢ — N) cos 6i' + sin(¢? — N) cos 6j' + sin 6k. 
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a andc are therefore 


— sin tj’ + costk 


) 
il 


c = costj'+sintk. 


Portion of the Reference Plane 





a 
as 


Y 


Figure 5.4: Angle Relations 


Since a and b are perpendicular, a- b = 0, therefore 
tan é = tani sin(¢ — 9). (5.1) 


Continuing 


b-c = |b| |c] cos(@ — 1/2) 


results in 


sin(@ —N) cosé cost +siné sint = sin#. 
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And using the relation from equation (5.1) results in 
siné = sin@ sint. (5.2) 
Finally, from c x b, the following relationship is obtained 
cos 6 = cos @sec (6 — N). (5.3) 
(5.1), (5.2), and (5.3) are the required angle relations. 


D. THE POTENTIAL 

In Chapter III, a simplified gravitational potential GM/r was used with the 
assumption that the attracting body was spherically symmetric. This simplified 
potential caused the satellite to move in conic orbits. As has been stated, the 
planets are not spherically symmetric but are bulged at the equator, flattened at 
the poles, and are generally asymmetric. An expanded expression for the gravi- 
tational potential will be developed in this section. The final expression for the 
potential is subject to the assumption made in Section B of this chapter that the 
planet about which the satellite revolves is approximately an oblate spheroid. 

Now, regardless of the nature of an attracting body, the potential V must 
satisfy one of the following differential equations. For regions within the attracting 
matter 


V’V =47pG where p is the density (5.4) 


which is Poisson’s equation. 


For regions outside attracting matter 
V’vV =0 (5.5) 
which is Laplace’s equation. 
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This discussion will be restricted to Laplace’s equation which, if V is written 


as a function of the spherical coordinates (r,0@,¢), may be written as 


18, Te av 1 av 
— —_— —_=—— SEE i LOE os = Q. 5.6 
r? Or (> M+ ates (co or) + ara om 


Any conservative force field F can be written as the gradient of a potential. 
Equations (5.5) and (5.6) state that div (grad V) is zero. The assumption that F 
is rotationally symmetric means that OV /0¢ is zero. Therefore a solution of (5.6) 


is sought that is a product of a function of r alone and a function of 6 alone: 
V (r,6) = fi(r) - fo(6). 

By multiplying (5.6) by r?, dividing by f,(r)f2(5), and rearranging terms, the 
following result is obtained. 

ld > df 1 O df 

—— | r°—] = ->———— — | cosh — ]}. 

fi dr dr focos 6 06 dé 
Since the left side is a function of r alone and the right side is a function of 6 


alone, both members are equal to a constant, say n(n + 1). f; must satisfy 


< (22) =n(n+1)f; 


which has the solution 
fr, = Ar” * 2ae,”. 
It is desired that this function be zero at r = oo, and that it be analytic and 
single-valued there. Therefore n is chosen to be an integer greater than zero, and 
B most equal zero. Therefore 
f= 


The equation for fe is 


d df, 
a5 (cos +n(n+1)coséf, =0 
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which may be rewritten as 


d d 
dain’ a Tm en ae (5-7) 





This is known as Legendre’s equation. Using the method of Frobenius, a series 


solution to Legrendre’s equation is found to be 
fo = CP,,(sin 6). 
Therefore the solution to (5.6) may finally be written as 
V(r,6) = Ar-” 'P,, (sin 6) 


where C is chosen as 1. 


There are many ways to consider the polynomials P,. Rodrique’s formula for 


P,,(sin 6) is 





ec 
P, (sin 6) = ———~—_(sin? 6 — 1)". 

(sind) = senld(sinayrom™ O—Y) 
so that 


1 
eo eee Sino, a — 5 (3sin’é —1), etc. 


The general solution of (5.6) is then written as 


ea) = s Anr~"~’P,(sin 6). (5.8) 


n=0 

The above mathematical derivation of an expression for the potential may 
be enhanced by a more physical formulation. Referring to Figure 5.5, let a unit 
mass m be placed at a point P which is a distance r from the center of mass of 
a bounded distribution of total mass M. Let dm be an element of the mass at a 


distance from O. Then the potential at P due to dm is 


—K*dm 


a —— 
(€2 +r? — 2récos 6)1/? 
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and the total potential is 


ge i (€2 + r? — 2récos6)'/2 (C2 


p 





Figure 5.5: Potential at an Exterior Point Due to an Irregular Mass 


The denominator of the right side of (5.9) may be expanded such that 
2g? ap 1a fey" . 
(é° + r° — 2récos @) ey = P,,(sin 6) 
a n=] 
where 6 = 17/2 — 6 


Therefore the following may be written 


_ 
y = K 





Af am tsinddm + A dm 


3 2 
Sf. ? «inne 


A description of the physical significance of each of the above integrals can be 
made. The first integral is the total mass, the second is the first moment about 
an axis through O perpendicular to OP and is zero when the origin is chosen as 
the center of mass. The third integral is the moment of inertia about the origin, 


and the last integral may be written as 
i (0? — f° cos* 5)dm = Ip — I 
M 


where Ip is the moment of inertia about the origin and J is the moment of inertia 


about the line OP. 
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The potential is then 


K’M_ K? 
A — = 573 (240 = 3). 


Now the first term is the potential due to a homogeneous solid sphere. The second 
term arises from the departure of the mass M from spherical shape. 


With the physical description established, (5.8) may be written as 


y=-<* a+ J, (2) (5 - Gin?) +... (5.10) 

Note that only Po and P,2 remain. 
An advantage of (5.10) is that it can immediately be written down once axial 
symmetry has been assumed, and any suitable experiment (such as an orbit of a 
satellite about the Earth) can be used to determine J, (or J3, J4, etc.; for genera! 


potential). 


E. THE EQUATIONS OF MOTION 


Referring again to Figure 5.1, the components of the velocity of a satellite in 
the coordinate system are: 


v= Ue ae as 


OT 


An expression for the kinetic energy (T) is 


GP \ >...» Conan ae on os ae" 
ae a 6 . 5.11 
oe (= eet (+) + r°cos ai | (5.11) 


The equations of motion may be written using Lagrange’s equation 
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da(rT-V) oa(T-V) _ 


= 9.12 


where g = a. T is the kinetic energy and V is the potential energy. 


Applying (5.12) to the expression for T (5.11) and V (5.10) results in 


d?r dé\° ,.(d¢\’ = av 
a -*(S) Be 5(3) ae a 
d{,d6\ .. d¢\* av 
ra a) sin 6 cos 6 ( $6) ae (5.14) 
d fie. 2,4b\ _ 
7 ( cos te) = 0 (5.15) 


The last of these equations gives an integral which can be used to eliminate 


time from the system, i.e., 


2 do 


oe cart e 
¥- COS O-—— = 71 COS a 
dt ° 


Let r = p/u then 
d h cos tou? sec?6 d 


Starting with equation (5.14) and applying (5.16) to the first term results in 


dt \" dt} dt p? do’ 
and continuing to apply (5.16) so that (5.14) finally becomes 


d ( S) d (= cos? to nee dé 





d* tan 6 We. 5 OV r*cos* 6 
——— + tan 6 = ——_- —__ 
dg? O06 h? cos? to 
or 
d* tan 6 sin 6 cos® 6 
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where 


3 R? 


p=h’/GM. 
It is desired to get the left side of (5.17) in terms of the independent variable 


6. To begin, the right side of the angle relationship (5.1), 
tan 6 = tan 1 sin(¢ —N), 


is substituted into (5.17). 
Then using (5.3) 


cos 6 = cos@sec(¢ —7N), 


(5.17) becomes 


tan 6 is an _ 2sec*1 cos 8 ,_2 v tan 2cos 6 i = Oe 
e cos 6 ¢'} ¢! cos 6 gp! 


aN? 
t & gto 
—2sec’1 tan 6 (= alas (| 


¢! sin 1 cos 2 g!? 


sin 6 cos’ 6 
= rs — a (2a) u), (5.18) 


where primes denote differentiation with respect to @. 

Later, ¢ will be eliminated from (5.18), but first the remaining equation of 
motion (5.13) will be rewritten in a form like (5.18). 

To rewrite equation (5.13) in a useful form, the independent variable is again 


changed by use of (5.16). Multiplying the expression by (¢')? and noting that 


(tan? 5)' = 2 tan 6 sec’ 5 
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results in: 


F 2.9, , (tand)”} | u(tan?6)' — u'd" 
u + uld*cos*é+ (ceca 5) | a. or 
a 
g'cos* 6 OV (5.19) 


h2u2 Or 





2 
Evaluating ae and multiplying (5.19) by con (4) results in: 


1 


ao, 
~[1 + J u2(1 — 3 sin? 6)] + u" 1- = ; 
(1 + seed —?'cos@sin @ tan i] 


COS 





i. i a 
COS“ to 


t Zs 4 fl 
u! cos? sect 6 E sin 6 cos 6(tan 6)! — q 


gi? 
t 6 12 O De 
— wu cos 1 sec? 6 + ee = ] , (5.20) 


Equations (5.18) and (5.20) are the equations of motion. Before proceeding 


to solve these equations, it is necessary to remove ¢. Combining equations (5.1) - 
(5.3) results in the following expression: 


tan(¢ —M) = cos 1 tan 8. 


Differentiating this expression with respect to @ results in 








nee COS 1 a ’ sin @ cos @ sin zt (5.21) 
~ cos? 6 cos? 6 
or 
1 cos’6 1 ; 
tae : yocos?é sy ; (5.22) 
@ COS 2 1 ?' sin @ cos 6 tan 1 


From (5.2), sin’ 6 = (sin 6 sin 1)? or 
cos’ 6 = 1 — (sin 6 sin 2)’. 
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Substitution of this expression plus (5.21) and (5.22) in (5.18) and (5.20) 
results in completely general equations of motion in terms of the dependent vari- 


ables: 1,,y; and the independent variable 6. 


ds e "3 “ 
6 sin t (%) sin” 0 ae sin 14% £! cos 6 sin? 6 


dé doz 
COS 2 COS t 
in 3 sin’? et 2 cost 6 sin 8 i 6 sin? t (29° cos? 6 sin 8 
COS 2 COS t 
re 3 sin i Ge cos? @ sin 6 i 2 sin’ 15 Ga? Sin g 
COS 2 COS 2 
2 di d?N 
, ee , eran 
"0 dg? COS 2 
3 sin >; (a0) sin 6 wdGane 
+ —————+—-_ + 6 cos sin: | — =] sin 6 
COS 1 dé 
2 e 
- 3 sin 1 (3) sin 0 ff sin? +23 40 sin 6 
COS 2 COS 1 
di dO d7i dO sin @ dQ?) 

3 co sin 9 — 2248" __ 2 sin : — sin 0 
an St oA q6 sin aaa Sia sin 
- ce 6 sin i? sin 6 3 sin’ 122 £0 cos? 6 

dé? COS 1 COS 1 

12 sin? 1222 cos? 6 dQ dn 
ats eased doa 3 cos 1 sin :— s 6 

COS 1 dé ge 
an 2 sin? 1222 cos 6 
eh — Cos, 6 2 
dG? COS 1 
di dP B229 cos 6 di 

eos — cos 0 ee cos 0 

er emeOS log do ae d0 


_ 2cos*i sin tJ u sin 8 (5.23) 


cos? i} 


and equation (5.20) is 
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d7u 3 sin? (#)" ie sin?(26) 
do? ia 4 cos?2 


3 -di dO d7u .-_ 2 ° 


cos? 2 


83 - di d?N du 
sin’: <3 & sin’ @ sin(26) 
2 cos? 


- 3 di dQ du ; : 
2 sin” 14,5, 9; cos @ sin @ sin(26) 


cos? 


sin’? (#)°u cos @sin @sin(2z) sin’: €) & cos? @ sin(2z) 


= : 
cos? 2 2 cos*2 
3 di dQ 2p: ‘di dQ du © 
sin” 15, 9,u cos’ 6sin(26) " 3810 (ae ae eine} 
cos? 1 cos? 2 
di d?u u dt a?7N N du 
sin 15 95 sin(20) i sin 15S & sin(2@) 
COS 2 2 cos? 
an d7u ae | 
3 sin‘ t (22) ao? on G _3 sin’? oa = au sin 49 
COs? 2 cos? 2 
ds dM du 2 di d7i du 
2 sin® ga & sin* eee Sor arr *9 
cos? 1 cos? 2 
- 3 -di d?0 du - 3 
2 SIN arr) C0500 ST aU 
cos? 1 
6 sin‘? (20 a cos 6 sin*@ 
cos? 2 
3: Pa iad } dD du 
2 SiN tage ap aCe a 39 
cos? 1 
2: 
2sin°2 Gr a cos 6 sin® 6 
cos? 1 
3 4 a dQ du 2 
9 sin a aes O° 6 sin? 6 
cos? 1 
- 2 di di du 2-di2 20 .:-.2 
2 Site. 7s an 60> eo rcine a sin“t7, u cos’ @ sin’ @ 
cos? 1 cos? 1 
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d7u oj2 2 ‘ 
6 sin? i (s0)’ gor in g _ 2 sin’: )d0 du cin? g 


‘dé do2 
cos? COS t 
in? 90 070 du 2:2 2 -d? du 
6 SIN ae asin 6 _ sin Pat as on 29 
cos? 1 COS 1 
8 «ds dQ du di dQ du 
2 sin”? eh oe og Sia “8, 2 sin 135 ae aa SiN 6g 
cos? 4 cos? t 
“i * 2 3 
Edt dl ed 2p Gt + de sin? g 
sin 1:— — —sin’ 9 — #200 dé 
d6 dé dé COS 2 
di d*t du . 4, sin ao a6 Sin’ 6 
— — —sin’ 9 — —“oe—_—_ 
d6 dé? dé COS 2 
3cos?tsin?:J usin? 6@ ~~ di? 
aa a sin? 6 
COS* to dé 
3 -d2i dO du SQ: 
SiN” 1353 ap gg COS 9 sin @ 
cos? 2 
8 sin?? @ a cos? @ sin 6 @ au cos® 6 sin 6 
cos? 1 cos? 2 
di du 4 sin’ iH oy cos? 6 sin @ 
— } —cos Se 
dé; dé Cos? 1 
‘di d?0 du : 
2 sin 17; 4,3 4, cos @ sin 0 
cos? 1 
dn \* du 
6 sin’ ce “aicOs 8 sin 8 
cos’ 1 
+ d24 dQ du 
2 sin 175%, 4, cos @ sin 0 
cos? 1 
d*i dO du 4 4 
sin 1—~ — — cos 6 sin 
dé? dé dé 
2 sin? 142 du cos 6 sin 6 _ sin ot & cos @ sin 6 
COS 1 COS 2 
2 sin?1 (st) 24 cos 6 sin 6 aaa 
dé) dé 
eS ee — eee Caen Sil 
cos? 1 d6/ dé 
— di df 2 sin iZu cos 6 sin 8 
4 sin 1— —u cos @ sin 6 — ———_—s—_ 
dé dé COS 1 


ol 


dn 
sin” 1% 22 24 cos ‘Oa 3 sin‘: (42 


a dé dé dé 


)" u cos‘ 6 


cos? 1 cos? 1 


sin? 74 20 44 cos? 9 


6 sin := 
i d6 dé dé a 


:dldu, 
+ 3 sin :— — —cos’@— 


d6 dé dé 


920 du _ 3a d@0du- d?M du 
4 446 dé? dé dé? d@ , dé? dé 


COS 1 cos? 1 COS 2 
dt dQ du di dQ du 


2 sin 1= 
_ 22 ae eee 


cos’ t d@ dé dé 
cos?1J u? ~—s cos? 2 


Cos? ty Cos? to 


3 sin? 


i—7u COS 


2 - 2 an 2 3 aa\? d?u 
- 3 sin’ (F) 1 CON rg Oe 3 (a) at 


dé 
COS 1 cos? 1 


-di dQ du 2 
6 do do cos’ 6 


cos? 1 cos? 2 
ds du 


£4 cos? 6 


dé dé 
COS 1 


dé d6? 


ds du 
sin 16 46 


COS 1 


(5.24) 


Equations (5.23) and (5.24) can be greatly simplified if only a first order 


approximation to the equations of motion is desired. The approximate equations 


to (5.23) and (5.24) are respectively: 


> 


t 
—2sinti—sin@ — sin 6 + sin :——cos @ — 2— cos 6 
dé dg? d@? dé 
2 cos*tsinitJu sin 6 
Cos? to 
d?u _ sin it TY sin (26) 
7 ai Aileen aaear mee 
dé COS 1 
2 ae: 
2 sin 20d sin?@ sin? 128 a sin?@ sin pad sin’ 6 
COS 1 cos? 2 COS 12 
3 cos?i sin?:Ju?sin?@ 2 sin’: oY -- cos 6 sin 6 
cos? COS to 
- d? t du 6 6 2 “dt Q ° 6 
sin tat ap eos sin sin 14 u cos 6 sin 
COS 1 cos 1 
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2 


(5.25) 





e - di du Z ° 2 °dQ 2 
7 3sint7,9,cos°@ 2 sin 17 u cos” @ 
COS 2 COS 1 
dQ d?u d? du > «di du 2¢ 2 
4 od0 a? , do? a6, SM tag 49 , COS tJu 
COs 1 cos 1 cos 1 COS? to 
cos? 1 
+ oe (5.26) 
COS” to 


Equations (5.25) and (5.26) will be used in the initial analysis; the general 


equations will be used for second order calculations. 


F. THE INITIAL CONDITIONS 

It is desired that the solution to the equations derived in the last section 
be an accurate, long term predictor of the satellite’s motion. In fact, as long as 
the oblateness perturbation remains the dominant disturbing force, the solution 
should be valid for close to 1000 revolutions. However, before the solution can be 
used for prediction, a set of initial conditions must be determined. 

The subject of orbit determination represents a separate discipline within 
celestial mechanics, and a discussion of the various techniques used to determine 
an orbit is outside the scope of this analysis. It is sufficient therefore, to state that 
the purpose of orbit determination is to find the orbital elements of a satellite from 
reduced observational data. A set of observations will determine the osculating 
elements at time tp. As was noted in Chapter 4, these elements will change, and 
at t; a new set of osculating elements may be calculated. For the purpose of this 
analysis any observed set may be designated as the osculating elements at to and 
thus prescribe the initial conditions. Stated mathematically, the initial conditions 


are at t = lo: 


> a(1 — e’) dr _ a(1 —e’)e sin(% — w) 
aie cos(@) — w) d6 (1+. cos(6) —w)? 


o3 


where a = ap, € = €9, w = wW, 


and t = to N= No y=6—-—w (5.27) 
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VI. THE PERTURBATION PROCEDURE 


A. PRELIMINARIES 

The perturbation method used in this analysis is a variation of the technique 
known as the Method of Strained Coordinates. The motivation for the use of this 
technique is the subject of this section. 

In perturbation theory, the quantities to be expanded can be functions of one 
or more variables besides the perturbation coordinate. An asymptotic expansion 


of f(@; J) in terms of the asymptotic sequence e¢,,(J) is, 


Ms 


f(0;J)~ Di an(en(J) as J+ (6.1) 


0 


m 


where @ is a scalar (or vector) variable independent of J and the coefficients a,, 


are functions of @ only. This expansion is said to be uniformly valid if, 


n—l 


do am(8)em(J) + Rn(9; J) 


m=0 


FOO lee Oe AP 


| 


f (4; J) 


For these uniformity conditions to hold, a,,(8)€,,(J) must be small compared to the 
preceding term dm_1(8)€m-1(J) for each m. Each term must be a small correction 
to the preceding term regardless of the value of @. 

Unfortunately, it is the rule rather than the exception that expansions like 
(6.1) are non-uniformly valid and break down in certain cases. A case of particular 
interest is the presence of secular terms such as 6" cos @ and 6" sin 6 which make 
Sm(8)/fm-1(6) unbounded as @ approaches infinity. 

To illustrate how secular terms arise in the solutions to differential equations, 


reference is made to one of the equations of motion derived in the preceding 


D9 


chapter. Equation (5.26) has the form, 
u"tu=14J f(u,u yu’, i,2",9',0"). 


When expanded, the right side will be linear combinations of trigonometric terms 
and constants. The presence of terms of the form cos(@), sin(@), cos(3y — 26), 
sin(2y — 6), etc., on the right side will produce secular terms. For example, the 


second order differential equation, 
il 
u +u=cos @ 


has the solution, 


8 sin 6+ cos 6 
a 


+ K,sin 6+ Ky cos @. 
Note that the presence of the term @ sin @ will result in unbounded solutions for 
u as § approaches infinity. Since in this representation wu is the reciprocal of 
the distance from the center of the planet to the satellite, the physical effect of 
the secular terms would be to produce in r periodic terms with large amplitude 
variations, a situation that certainly does not correspond to physical reality. 

A technique for dealing with terms that produce secular terms is to eliminate 
them from the right side of the differential equation. The Method of Strained 


Coordinates is a perturbation technique designed for removing secular terms. To 


illustrate, it is recalled that the solution to the two body problem is 
1/r =u=1+e cos(6 —w), where p= 1. 


In the Method of Strained Coordinates, the @ coordinate is strained by introducing 
a new variable y = f (@) = @—w+Jy,6+... As will be shown in the next section, 
choosing the correct value for y,; will insure secular terms do not arise in the first 


order solution for u. 
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While the above method removes secular terms to first order, it is not ade- 
quate for dealing with secular terms that arise to second order in the equations of 
motion. That it is important to remove second order secular terms is shown by 


the following equation: 


f(r) = 1+€ecos(6—w) + J(cos 6+...) 
+ J?(6 cos(6 —w) +...) + J*(6? cos(@ —w) +...) 4+... 


Note that, although the terms through order J are bounded, as 6 — oo, the J’ 
and J* terms grow without bound and dominate the right side. In this present 
analysis, 6 has an upper bound of (27)10°; however, the effect of secular terms 
remains since J76, J*6’, etc., are all of order J as 6 > (27)10°. An infinite series 
would have to be retained. 

In this present work, three additional techniques had to be devised to deal 
with secular terms to second order. These techniques will be discussed in the next 
section. The perturbation method therefore is not strictly the Method of Strained 
Coordinates, but a variation of that method. 


The basic steps in the perturbation procedure are as follows: 


1. The dependent variables and independent variables are expanded as func- 
tions of a small parameter (J). 


2. The variables are then substituted into the equations of motion, and the 


equations are solved consecutively. Each solution yields a more exact ex- 
pression for the appropriate variable. 


The process is carried out through second order to demonstrate that all secular 
terms may be eliminated and that the solutions are bounded. The following sec- 
tion highlights the calculations involved in the process and shows the first order 


equations and solutions. The second order expressions are long, and their display 
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in this context would not contribute to the analysis since only a few specific terms 
are relevant. The complete expressions are contained in Appendix B; they will be 


referred to during the course of the analysis. 


C. SOLVING THE EQUATIONS 
1. First Order Approximation for : and /2 
The equations to be solved are (5.25) and (5.26). Since the right side 
of these equations are analytic functions of J, it is reasonable to expect that the 
solution to u will be arbitrarily close to the two-body solution, 1 + e cos(@ — w), 
when J is sufficiently close to 0. Likewise 1, 9, and y will be arbitrarily close to 


19, No, and 6 —w, respectively, when J is close to 0. This assumption amounts to 


letting 
u = l+ecosy+Ju,+J*uz+... (5.2) 
y = 0O-wt+Jy4+J*y2 4+... (5-3) 
1 = igo t Ji, +J*t24+... (5.4) 
9 = 29+J0,+ 5724+... (5.5) 


The first step in the solution of the equations of motion is to substitute 


(6.2) - (6.5) into equation (5.25) and equate the coefficients of J. The result is 
(—20), sin 6+} cos 6) sin tp — 27, cos 6 — 1 sin 8 
= 2 os tgsin tgsin 6(1+€ cos y). (5.6) 


Sin t and cos 2 have been replaced in the above equations by their approximations: 
sin tg and coS to. These are valid approximations since M”, N’', 7", and 7’ are all of 


order J. 
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Equation (6.6) is a linear differential equation with two unknown func- 
tions. Terms of the form cos 9, sin 6, etc., on the right side of this equation, and 
the more general equation (5.23), will produce secular terms int and N. It should 
be recalled that these same terms cause secular terms in u. Note that (6.6) con- 
tains a sin @ terms on the right side. There would be a need to eliminate this 
term were it not for the conditions placed on NM and 1 by the definition of the 
reference plane. Specifically, 92, which governs the rotation of the reference plane, 
must be expressible as an unbounded (secular) term plus bounded periodic terms, 
and : must be bounded. The fact that M must contain a secular term negates the 


requirement to eliminate the sin @ term. Therefore, a solution of the form 


9, == a4 + asin y 


vy Bi cos y (a7) 


is assumed. 


By substituting (6.7) into (6.6) and equating coefficients, the following 


solution for 2 and 2 is obtained: 


—2/3e cos igsin tgcos y 


vy 
1, = —@ cos to —4/3e cos tosin y (5.8) 


2. Second Order Secular Terms in: and 2 
Equation (6.8) satisfies (6.6) to order J. The next step in the process 
is to substitute (6.2) - (6.5) and (6.8) into (5.26) and solve for u,;. However, 
proceeding with (6.8) in its present form will lead to secular terms in second 
order. A brief paragraph will explain why secular terms are anticipated. 
Success in using perturbation methods requires a certain a prtort knowl- 


edge about the nature of the particular problem one is trying to solve. There is a 
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certain amount of trial and error involved in the process. For instance, the success 
of the Method of Strained Coordinates depends on the knowledge that a secular 
term will arise in the first order solution to (5.26) and that a mechanism (the 
strained coordinate) can be put in place one step ahead in order to eliminate the 
secular term. The second order secular terms were found by this same trial and 
error process. 

Returning now to equation (6.8), it was discovered in this analysis that 


the problem terms 
(—45s° + 28s) 


7e7c sin(2y — 6.9 
Fi J*e*c sin(2y — 8) (6.9) 


and 

9 72.3.2, -: 

ao s“e’c sin(2y — 36) 

(where s=sin 19, ¢ = COs 129). 
appear on the right side of (5.23) when it is expressed to order J?. The appearance 
of these terms that give rise to secular terms was unexpected. They were not 
reported by Brenner because the authors assumed a small eccentricity and dropped 
J*e* and higher terms. 
There were several failed attempts in dealing with these new terms be- 

fore an effective measure was discovered. First, an investigation was made into 
the effect of retaining the secular terms. 


The secular terms produced by (6.9) are: 


7 15s° — 14 
—— ues’ — Mea 6 sin(2y — 26) 
24 
15s? —7 
a Cos 6 cos(2y — 26). 


As was demonstrated earlier. when 6 — (27)10° these terms are of order J and 


must be retained in the first order solution. As the solution progresses, these terms 
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will continue to produce secular terms with coefficients J°*6?, J405, ..., J"@"-). A 
convergent series representation for these terms was not found. 

The remaining alternative was to alter the form of (6.7) to eliminate 
secular terms. If terms with arbitrary coefficients can be found which when added 
to (6.7) produce terms of identical harmonics to the same order as (6.9), the co- 
efficients may be chosen so that all terms with these harmonics are eliminated. 
Essentially, the only terms that may be added to (6.7) are terms which satisfy the 
original differential equation, equation (5.23). Therefore, one may add homoge- 
neous solutions of the differential equation or arbitrary constants. 

Adding an arbitrary constant to (6.7) has the following effect on higher 
order terms. From equations (5.23) and (5.24), it would appear that only the 
derivatives of : and | enter into the higher order calculations, and therefore arbi- 
trary constants would be eliminated. This is true for N, but not for 1. To explain, 
it is recalled that the approximations sin 1 = sin 19 and cos t = COs to were valid 
for the first order approximation of 1, and ,. This approximation was valid since 
all terms were of order J. However, in the calculation for u,; there is a term 
cos’ 1/cos? i that is of order 1 (Eq. (5.26)), therefore a better approximation for 
1 is required. Using (6.7) and the Taylor series expansion for cos : (keeping only 


terms of order J) results in 


cos(tp + f(J)) = cos tp — sin tof (J) 


COS 2 


COS tg — SiN to(1y) 

Adding an arbitrary constant K, to (6.8) results in 
COS 1 = COS tg + sin 19(2/3 Je cos tgsin tg cos y + Kg) 
(A similar expression is required for sin 1) 
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Therefore, a constant added to (6.8) will alter the form of subsequent 
terms, but not in the form required to eliminate secular terms in 22 and M2. It will 
be shown that a constant like K, will produce a term of the form K,J sin’ to in wu; 
and terms of the form K,J’e sin? igcos tgcos y in tz, and K,J%e sin? t)cos tgsin y 
in M2. Although they cannot be used to eliminate secular terms, the constants 
will be needed to satisfy the initial condition, therefore it is essential that they do 
not produce irremovable secular terms to higher order. 

The next alternative for the elimination of secular terms is to add a 
solution to the homogeneous equation of (5.23). While it will be shown that this 
technique is successful in eliminating secular terms in ug, it did not succeed in 
removing them in tz and M2. It failed because the homogeneous terms produced 
new secular terms to higher order. Many various combinations of homogeneous 
solutions were added to (6.7), but all attempts along this line only complicated 
the problem. 

An answer to the question of how to eliminate the secular terms was 
suggested in a report by Weisfield [Ref. 31] on polar orbits. When faced with the 
problem of eliminating secular terms from his equation for A*, Weisfield added a 
term like cos(2y — 26). Now to the particular order to which one is working this 


term acts like a constant in the derivative, e.g., let, 
y=6+4+J86 
then 


© (cos(2u — 20)) =sin(2y — 26)(2 — 2(1+ J)) = O(J). 
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The attempt was then made to apply Weisfield’s technique to this more 


general problem. (6.7) would then become 


ty —2/3e cos igsin tocos y + K,cos(2y — 20) + Kz 
N; = 6 cos to — 4/3ecos tosin y+ Kgsin(2y — 26) + Kg. (6.10) 
where: K, =—Kye?cs Kz = K,cos(2w) + K; 
Ks =—Kse*c Kg = —Kgsin(2w) + Ky 
(It is noted that (6.10) still satisfies (6.6) to order J) 
The added harmonic terms produce identical harmonics as (6.9) to sec- 


ond order and no other new secular terms. Thus, they allow for the elimination 


of these problem terms. 
3. First Order Approximation for u 


With valid expressions for 1; and 1, established, the procedure may now 
be continued with the calculation of u,. Substitution of (6.2) - (6.5) and (6.10) 


into (5.26) and letting sin i9 = s and cos 19 = ¢ results in the following equation 


d’ 5e?s? cos(2y + 20 
Sen te = 2e cos yy; +o 
e*s* cos(2y — 26 
+ 2e?K,s’ cos(2y — 26) + eaeeesiey — 20) 
Ses’ cos(y +20)  11e’s? cos(2y) A 5e’ cos(2y) 
3 4 2 
3e’s? cos(26) if s’ cos(26) 
4 2 
Iiewoen 5S” . 7e° 


ie eae le GOmel 
12 ee ey ( ) 


5es’cos y + 4e cos y+ 





— 2e*?K,s* cos(w) — 2K2s* — 


In the above equation, the cos y terms will produce secular terms in u,. A choice 


of y; = (5s? — 4)/2 will eliminate that possibility. y becomes 
D2 2 
y=0-w 40 (5s -2)+J yf. 
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Equation (6.11) becomes 


d*u, 5e’s*(2y + 26) oa 
We? +u, = re + 2e°K,s° cos(2y — 28) 
e*s? cos(2y —26)  5es*cos(y+20) 11le?s? cos(2y) 
en + _—$._ 
8 3 4 
5e*cos(2y)  3e?s*cos(26) —s* cos(26) 
2 a oi 


17e*s* 58% = Te? 
et ee. (6a) 


2e*K,s* cos(2w) — 2K.s* — 
€ 4418S cos( Ww) 2S 12 9 6 





To solve (6.12), a solution of the form 


ur = a, cos(2y + 26) + a. cos(2y — 20) + as cos(y + 20) + agcos(2y) 


as cos(20) + ao (6.13) 


is assumed. 
Equation (6.13) is substituted into (6.12) and the coefficients of like harmonics are 
equated. The coefficients are: 


—e*s? 


24 





Aa = 


1 
= (2K, + = | es? 








—5es? 
a = 
} 24 
alls? 85, \. 5 
as = apes iG 
12 6 
«> -(G) 
a Vaaee 
lige, 5 
Qo = —2Ky,e’s’ cos(2w) — 2K2s’ — ae s? 
12 2 
peg 
wae 
6 
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Equating (6.13) for u; results in: 





e*s* cos(2y + 20 
u = —£ cos(y + 4) + 2e7Ks’ cos(2y — 26) 
e’s* cos(2y— 26) Ses? cos(y + 26) ‘a 1le*s? cos(2y) 
8 24 12 
Se? cos(2y) e?s?cos(26) —_s? cos(26) 5 
- oa 2e7 Ks’ cos(2w) 
Wes’ oS" tee 
a ee 2 esp | 


4. Second Order Secular Terms in u 
The above expression for u; is not complete. Looking one step ahead 
it is known (by trial and error) that the following problem terms arise in the 


differential equation for up. 


d*uz i 75e°s* — ii 75e*s* — 200e*s* + 136e° e° (3y — 26) 
— = ———— } cos(3y — 
dg2 ~ 15(5s? — 4) : 
375e°s4 480e* — 40 136 ‘ 
| SE ek Be) ee FD ng — 20 
240 15(5s? — 4) 
(225e°s* — 225e°s? + 2e°) cos(2w)  2e* cos(2w) 
—- | 2ey. & +. —_____. 
60 75s? — 60 


(45e° — 550e)s* + (36° + 488e)s° — 56e 


- Bea ee eee | cos(y) (6.15) 
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(Only the problem terms are displayed. The complete expression is equation (B.5), 
App. B.) 

The harmonics cos(3y — 26), cos(y — 26), and cos y all cause secular 
terms in uz, and therefore must be eliminated. In addition there is another more 
troubling problem with (6.15). That is, certain terms have in their denominator 


(5sin??) — 4), and if 1) = +sin7? \/4/5, then the denominators are zero. This in- 
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clination is well known in Astrodynamics, and has been named the critical inclina- 
tion, or more appropriately the critical inclinations since there are two: t9 = 63°26' 
and 117°34'. The problem of the critical inclination will be dealt with later in the | 
analysis. In addition, there is qualitative discussion of the critical inclination in 
Appendix C. 

The task is to eliminate the three terms which give rise to secular terms 
in u2 from the right side of (6.15). By inspection it is seen that a proper choice of 
y2 will eliminate secular terms produced by cos y. It was then discovered in this 
analysis that the addition of a term of the form cos(y— 26) —cos(y+2w) eliminated 
the cos(y — 26) problem term. This term may be added since it is a solution to 
the homogeneous equation. There remained one term to eliminate; cos(3y — 29). 
This term can be eliminated by adding a term with the harmonic sin(2y — 26) to 
y,. In addition, a constant was added to y; to satisfy the initial conditions. The 


complete first order expression for y is now 
5 See : 
y = pw I (8 (Fs - 2) + Resin(2y — 26) + Ke) 
+ J? y08 . (6.16) 


where Ks = —JKsand Kg = —JKssin(2w) + Ke. 
By adding — K7(cos(y—20) +cos(y+2w)) to eliminate a secular term and Kg cos y 
+ Kgsin y to satisfy initial conditions, u,; becomes 


e’s? cos(2y + 26) 


uu; = 3, + 2e*K 18" cos(2y — 26) 
*s* cos(2y — 20) 5es? 20 
ee cos(2y ) _ 5es* cos(y + 29) wre cos( 28) 
8 24 
1le*s? cos(2 Se’ cos(2 
+ eKy7cos(y + 2w) + a ; sd + Kgsin(y) 
a5 26 s 260 
+ Kegcos(y) — <s cosi?* — = sos(20) — 2e? Kis’ cos(2w) 
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17e*s* 5s? = Te? 
6 A ill (6.17) 


a 2K. s? a 
12 2 6 





5. Second Order Solution for : and 2 


With all terms in place to deal with secular terms the calculations are 
continued by substituting (6.2) - (6.5), (6.10), (6.16), and (6.17) into (5.23) to solve 
for MN, and t2, and to evaluate the constants K, and Ks. The result is contained 


in Appendix B, i.e., 


(—2N, sin 6+ 12 cos 6) sin to — 21, cos 6 — 13 sin 6 
= Right side (B.1), (Appendix B). 


An inspection of (B.1) reveals that the coefficients of sin(2y — @) and sin(2y — 36) 


form respectively the following simultaneous equations: 


(5Ks+10Ki—45)e*s°c  (4K3+4K,—7)e*sc _ . 
24 6 - 
0 2.3 2 
3 7 OKs e’s°c + (4K,—4K3)e*sc =0. 
Solving these equations results in 
= = 15s? — 14 (6.18) 
' “ 24(5s? — 4) ] 
75s* — 120s* + 56 
k; = —— . 6.19 
: 24(5s? — 4)? es) 


Substitution of these values for K, and K3 into (B.1) gets rid of all sin(y — 26) 
and sin(2y — 30) terms. With an assurance that no secular terms will arise, the 
equation for MN, and 7 can be solved. 

Before progressing, it is noted that (6.18) and (6.19) contain tne same 
problem denominator that was observed in (6.15). In fact, (6.18) and (6.19) are 


the first occurrence of the critical inclination term in this analysis. The term then 
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continues to manifest itself in all higher order analyses. The reason for the presence 
of the term is simple to explain. The constants K, and Ks were multiplied by the 
derivative of y; (4 = (Es? — 2)) during the analysis. The derivative then shows 
up in the denominator of these constants when they are evaluated. It will be 
important to show later in the results that despite the occurrence of an apparent 
singularity, the final solution is uniformly valid for all inclinations, including the 
critical. 

With secular terms removed the following equation can be solved for N2 
and to. 


(—2N) sin 6+ 1} cos 6) sin tp — 21, cos @ — t) sin 0 (6.20) 
= Right side (B.2), (Appendix B.). 


As was done in (6.6), a solution which includes the harmonics of the right side of 
(B.2) is assumed for NM; and t2. This solution is substituted into (B.2), and the 
coefficients of the various harmonics are equated. Once again the conditions are 
that 2, must be expressed as a secular term and bounded periodic terms, while 


12 must be bounded. The solution is 
NM, = Right side (B.3) 


Right side (B.4). (6.21) 


1 
6. Eliminating the Final Secular Terms 
To complete the solution the remaining constants must be found. Three 
of the remaining six constants, y., Ks and K7, are obtained from the differential 
equation for u,. It is recalled that these constants are used to eliminate secular 


terms in uz. There is no need to solve the equation for uz since the constants may 


be evaluated from the right side of the differential equation. Therefore, (6.2) - 
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(6.5), (6.10), (6.16), (6.17), and (6.21) are substituted into (5.24). The resulting 
equation is (B.5). 

As in the previous solution for NM, and t2, the secular terms in uw» are 
eliminated by setting the coefficients of the problem terms, i.e., (6.17), equal to 


zero and solving for the constants. The results are: 


75e7s° — 260e?s4 + 296e7s? — 112e? 





K; } 
192 (Ss? — 2) 
15e’s4 — (14e? — 2)s? 
a 
48 (Es? _ 2) 
e* cos(2w) 
2 oa — ape ae Nari (6.22) 
30 (Es? — 2) 
where 
26es? = 2.4 
Y2 = —10es*cos(w — 4) + 26es" cos(w — 4) _ 15e*s* cos(2w) 
3 8 
* 15e’s? cos(2w) e* cos(2w) if 15e’s4 27584 3e7s? 61s? 
8 60 32 48 g 12 
re 
lez 


With the above constants evaluated, it is assured that the second order expression 
is free of secular terms. 
7. The Initial Conditions for y, 1, 8, 9 
The task now is to evaluate the remaining constants by establishing the 
initial conditions. At t = fo, it is required that the velocity vector of the satellite in 
the reference plane be tangent to the corresponding two body ellipse determined 
by the satellite. In addition to the initial conditions established for r and a re 


(5.27), it is recalled that at t = to 


y= 9)-w, — a. and Ge= (1 
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From (6.16) and (6.22): 
° 
a p—w4 J (6 (25-2) 
200 2.4 296e?s? — 112e? 
(75e*s° — 260e*s* + a € ) (sin (2y ~ 26) sin(2u))) 
192 (Ss? — 2) 
e?J7@ cos(2w) 


4 (Es? — 2) 


+ J76 Go. (6.23) 


Choose Kg such that at t = to, y = 0) — w; therefore, 
= —J (52-2). 029 


To obtain the initial condition for t, use (6.10) and (6.18) so that 


(15s? — 14)(cos(2y — 26) — cos(2w)) 


1 = 1%) — 2/3Jecs cos y — Jees 24(55? — 4) 


+ K,. (6:25) 


It is evident that the Je? terms are eliminated when y = 6 —w and 6 = @). The 
result 
Ky, = 2/3Jecs cos(@ — w) (6.26) 
gives the desired 1 = to at t = to. 
For O, (6.10) and (6.19) are required. In addition all secular terms from 


M., equation (B.3), App. B, are needed since all these terms are of order J. 


Therefore, 


ce? J (75st — 102s’ + 56)(sin(2y — 26) + sin(2w)) 
24(5s? — 4)? 

4ceJ siny  ce?J* cos(2w)6 a 5ce?J*s? cos(2w) 6 

3 15s? — 12 8 

ce’ J? cos(2w) 6 Seems 6  5odas-6 
= ————— + 5c J? K, 876 — —__ 

12 isiaaitiis 24 (OS 
Ce Jiggmeecd 0 
6 Z 


9 = No-Jbe+ 











ame’ (6.27) 
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and 
4 
Ky, = Jc + yee sin(Oy — w). (6.28) 


8. §@ As a Function of Time. (Initial Conditions for r, 

To complete the evaluation of the remaining constants in u, it is nec- 
essary to give @ as a function of time. The expression for the perturbed “ then 
can be related to the well-known two-body formula, and from this relationship, an 
expression for r and a may be derived. The task will be to choose Kg and Ky 50 
the conditions as given by (5.27) are satisfied, i.e., 

ea = : a(1 — e?) =“ _ a(1 —e’)e sin(9 — w) 

+e cos(8> — w) dé (1 + e cos(8> — w))? 

Proceeding in this manner, the formula for a can be obtained from the 

relation 


dt dpdt_ d®¢ r? cos? 6 


dd d0dd dO hos io 


Direct substitution of the derived expression for a equation (5.21), results in 


dt ‘a nae _e’ cos(2y — 28) — e’ cos(2w) " e*s* cos(2y — 26) 
d0 h 15(5s? — 4) g 
e?cos(2y — 26)  es*?cos(y+20) es*cos(y — 26) 
60 6 Z 
in 4es’cosy 4ecosy  s*cos(20) 2es*? cos(w — 9) 
3 3 Z 3 
e*s*cos(2w)  e?cos(2w) — s? 
= + ——__— +—-]l 6.29 
8 i 60 i ( ) 
Peat — to 
2 


a 


The above condition is an integral from the two-body problem. If (6.29) is evalu- 
ated at 6 = 4, then # + may be expressed as 


r? 
(1+ JK) => 


sills 


where 


es?cos(w +6)  2es*cos(w— 4)  4e cos(w — 4) 
a -- Ne 
i 3 3 
es? cos(w — 38)  s*cos(28) 


cree 
6 2 BS 


1 


From this expression, a formula for h results: 


so that from (6.29) 


dtr’ 
Ta zl + J(-Eio + ---)) 
Using 
h? h? 


Daa e?)(1 + 20 Kyo) 


a new formula for r, which includes the first order perturbation effects, may be 


written as 


a(1—e’) 


p 
a 1+ecos y+ J(—2Kjo + u;) 
). 


(6.30) 
where u, is (6.17 


The formula for x is 


dr a(1—e*)esin y 


dd — (1+ e cos y)? 


Ja(e?—1) |[e*s*sin(2y + 26) ” 5es* sin(y + 26) 


(1+ ecos y)? 6 8 — eK7sin(y — 26) 
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é 1] 2 62 G} 22° 2.° 
= eK, sin(y + 2w) — lle*s* sin(2y) as Se" sin(2y) Es es: sin(y) 


Kgsin y+ 2e sin y+ Kogcos y+ 


6 3 2 
e*s* sin(26) - s* sin(2@) 


: ; (6.31) 


Next, evaluate (6.30) and (6.31) at 8 = 0. Keep only the terms to order J. The 


result is two simultaneous equations that can be solved for Kg and Ko: 


Kg 


Kg 


_e?s’ cos(3w — 0) | 11e?s? cos(3w — 360) 


16 24 
5e? cos(3w — 360) _ e?s? cos(3w — 54) HY 31les? cos(2w — 26 ) 
4 16 tz 
Te cos(2w — 269) _ 3es*cos(2w — 48) — 11le*s* cos(w + 8) 
3 8 16 
s*cos(w +69)  1le*s*cos(w — 69) | 7s? cos(w — 8) 
4 8 Z 
31e* cos(w — 17e’?s* cos(w — 34 
—- — 3cos(w — 4) — <a 
7s* cos(w — 39) es’ cos(2w) Ses’ cos(28) | 18es* Te (6.32) 
2 2 4 12 3 


e*s*sin(3w— 6) 11e’s* sin(3w — 36 ) Fi 5e?(3w — 38) 


16 24 12 
e?s*?sin(3w — 56.)  3les*sin(2w—26)) Te sin(2w — 24) 
16 12 3 
Ses’ sin(2w — 44) _ 7e*s?sin(w +o) | s?sin(w + 4) 

8 16 4 
67e*s*sin(w—O) 7s?sin(w—O)  29e?sin(w — 4) 

24 ‘ 2 vi 12 

: ae 
3 sin(w — 6) + Liets: nae — 34) és 7s sal — 36) 
es’sin(2w)  3es’ sin(24) | (6.33) 
Zi é 


fe 


Substitution of the values of Ag and Kg in the equation for u insures 


that the remaining initial conditions for r and a are satisfied. 


C. THE RESULTS 
1. Preliminaries 
It was noted that throughout the calculations in this chapter the term 
(Bein 0-2) appeared in the denominator of terms in tu, 12, M2, and uz. It would 
appear that the solution is not valid when tp) = +sin7’ /4/5. And if the solution 
is not valid near these “critical inclinations”, then one should question the validity 
of the underlying perturbation process. 
The purpose of this last section is to display the final solutions for y, 1, 
M, and u, and to show that each of the apparent singularities can be eliminated 
in the limit as as 19 approaches +sin™? /4/5. It is a remarkable aspect of this 
analysis that every term that appears to cause a solution to blow up is exactly 
canceled by a corresponding term that is “hidden” in the solution. 
2. First Order Solution for y 


By use of trigonometric identities, (6.23) may be rewritten as 


r— 
| 


5 
6—w+J0(=s-2) 


758° — 260s* + 296s? — 112 5 5 
— pS OS MN) (2. — Jé@ (<5? - 2) ) sin (Je (25? —2 


96 (£5? — 2)" 2 


e7J*6 cos(2w) 


o 
= + - S00 (25? — 2) ro) Oi iO Glial eta -)- (6. 


30 (§s? — 2) 


Equation (6.34) is the complete first order solution for y. If the initial 


inclination is such that |3s? — 2] < 107%, then (6.34) should be replaced by its 
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34) 


limit: 


e? J* sin(2w)6? is 8eJ* cos(w — 40) 8 m 17e’? J* cos(2w) 








= @-w+ 
v 15 15 60 
eJ7@ 2/76 5 
0 (J?, 1% 3*6° (=s? -2) we 
+ i + : + , J 16 5s Ze 


3. First Order Solution for : 


As was done with the expression for y, 1 is rewritten as: 
, 2 
1 = to+ Jes (F¢lcos(6o —w) —cos y| 


e?(15s” — 14) sin |2w — J0 (8s? — 2)| sin [70 (8s? — 2)| 
24 (Ss? — 2] 


+ O(J?+ 3°64 S10 +...). (6.35) 


Again, as in y, should 19 be near the critical inclinations, [257 =2| =< 16~, 


then (6.35) should be replaced by: 


J7*9e? sin 2w 


. . Ase 
1 = to+ 75 [c0s(8o — w) — cos(8 — w)] + 30 


5 
lO) (J?, J°0, 3°60 (= sin® i 2) +...) 


4. First Order Solution for 


Rewriting (6.27) for results in 


Q ee 


75s4 — 120s? + 56 
a 
48 (3s? — 2) 


ox (200 82-2) an (24 ($2) 


4ceJ siny  ce’J* cos(2w)é Ka 5ce*J*s* cos(2w)6 
3 15s? —a2 8 
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ce? J? cos(2w)6 See? J7s76 = 5eJ* 576 
tO fe ee 
12 OCS ae 2m.!~CS 

ce*J*6 =e 78 
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+ Ky+ O(J’, J°0, 5467 +...) (6.36) 
where 
2 
Kz, = 3Jecs cos(6> — w) 


4 
Kg = Jc0o+ vce sin(@) — w). 


If |2s? — 2| < 10-*, then 

4ceJ sin(@—w)  ce*J*sin(2w)6? A 5ce?J? cos(2w)é 
3 6 rr 

ce?J70 = 11078 
uM 6 


— 
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Q9 — Jbe — 








ae 4cJ*K26 — + K¢ 


O G ogee 6: (2s = 2) ce ...) 
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where 


Ky 


= Je cos(8) — w) 


4 . 
\/1/5 J0o+ 3V1/5 Je sin(9) —w). 


5. First Order Solution for u 


K4 


The complete first order solution for u is: 


e? J s* cos(2y — 26) 


1 eicosys a) as + 2e7J Kys’ cos(2y — 26) 
a e’ Js’ cos(2y — 26)  5eJs* cos(y + 28) as lle*Js*cos(2y) _ 2K; coaly aan 
8 24 2 
5e7J cos(2y) e*Js? cos(28) Js? cos(26) 


ar ed Hy cosy 2120) — a eee 6 


— JKogsin(w — 0) + JKgcos(w — 69) — 2e?J Kys’ cos(2w) — 2J Kys’ 


lie Js oo tee 
ee 
ee 2 











+ J + O(J’) (6.37) 
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where 











15s? — 
i —_ : 14 
48 (8s? — 2) 
Ky = _ 2ecs cos(@ — w) 
y, = p-w+J0(- —§ *—2) 
75s° — 260s‘ + 296s? — 112 
— et ee + ee 1?) s (2w — Ja (5 —S *—2)) sin (v0 (2 S *—2)) 
48 (8s? 2 2 
— e'J*@ cos(2w) + JA (35 — 2) + J709, (see equation 6.22 for 9) 
pe ptcd — 2 9)\¢2 
fe = 15e*s* — (14e° — 2)s 
48 (3 25? — 2) 
Kg -—— Equation (6.32) 
Ky —— Equation (6.33) 
Applying the condition |2s? — 2| < 107°, the expression for u is 
P _ e?J cos(40—2w)  e? J cos(2w) — eJ_cos(38 — w) 
— 30 10 6 
*J cos(26 — 2 
= nim + JKgsin(@) — w) + J Kg cos(8 — w) 
(<= sin(2w)6*  8eJ? cos(w — 60)@ 
+ e cos | ———————_ + 
15 15 
1 17e? J? cos(2w) 6 ” eee 19 
— —W 
60 60 5 
e(8e* — 8)J*A sin(0+w) e?J cos(20) — 2J cos(28) 
20 5 15 
272¢: 
det J Sie ial alaailale Pe (262 2) +. » 
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Where Kg and Ko are now: 


e? cos(3w — 8) 


Kg 
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e? cos(3w — 56o) 


20 20 


e* cos(3w — 3p) 


4e cos(2w — 26)  3e cos(2w — 469) lle? cos(w + 4) 


mem mS eee 


15 10 20 
cos(w +o)  89e’cos(w — Go) _— cos(w — Ao) 
5 60 5 
_ 17e’ cos(w — 38) _ 7 cos(w — 38) _ 2e cos(2w) _ e c0s(260) 
60 15 5 
22e 
—- —+40(J’ 
i 


e’sin(3w — 4) ii e’ sin(3w — 36 ) ad e? sin(3w — 56) 


a 20 20 20 
a 4e sin(2w — 26) . 8e sin(2w — 46) Te’ sin(w + 4) 
15 10 20 
A sin(w +69)  lle*sin(w—@)  sin(w — 4) i lle’ sin(w — 36) 
o 60 o 60 
ie 7 =— 36) : 2e aS _ 8e 4 + 0(J?) 


The results obtained thus far may be used to predict the orbit of a satel- 
lite for up to 1000 revolutions when a valid set of initial conditions are provided. 
The initial conditions will usually be given as the initial displacement vector r 
and the initial velocity v, whereas the solution obtained is in terms of the orbital 
elements. It is important to show that the coordinates and velocity components 
can be easily recovered from the orbital elements. Indeed, Reference [4] contains 
the criticism that Brenner and Latta did not show how r and v would be ob- 
tained, and that “the derived elements are not such that the velocity components 
are readily obtainable.” (Arsenault et al., Ref. 4: Vol. 2, p. 5). 

Again referring back to Figure 5.4 and equations (5.1) - (5.3), the po- 
sition of the satellite in the coordinate system may be derived from its direction 


cosines: 


r = rcosécos ¢i+rcos ésin¢dj+rsin 6k 
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where r is equation (6.30) and i,j,k, are unit vectors. 


Measuring from the line of nodes results in 


r = rcos 6 cos(¢—MN)i +r cosé sin(¢ —N)j' +r sin 6k 
where i’ = cosNi+sin Nj 
j = sin Ni+cos Nj. 


Using the angle relationships from equation (5.1) - (5.3), the above expression is 


rewritten as 
r = r(cos@ cos 2 —cos? sin 6 sin N)i 
+ r(cos @sin M+ cos? sin@ cos M)j+r sini sin 6k. (6.38) 


The velocity is found by differentiating (6.38) with respect to @ and noting 


dr _ dr dé 
dt d6dt~ 
The result is 


_ fdr, . dro, | dr3, \ dé 
v= (G+ S54 Bu) dt 


where 


d d df 
= —cos2 sin n= sin 6 — cos 1cos Nat sin @ 


di 
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+ cos 1-7, cos 8 sin QT cos 
— cosisin N2rcos @ 
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76 r ( cosi sin’ 70 sin 6 —sin 9 sin 


79 


. dt ; dQ 
— sin 77 cos 21 sin 6+ cos a cos 6 


+ cos:cos 11 cos 6) 


ae = (cos :cos M sin 8+ sin N cos 6) 


drs rsint cos 6+ sin ¢ sin gat +rcos1sin gu 
——= ‘= 7 2 $l os 
dé dé dé 
The requirement now is to find expressions for 4, 4, and %; the expression for 
* is equation (6.31). 
Differentiating the expression for :, equation (6.35), and the equation 


for N, equation (6.35), the following equations for a and a are obtained: 


< = 92) €cs simpy 
dQ} 4 . 
Tiss —Jc— gJce sin y. (6.39) 


The expression for é is equation (6.29). By expanding this expression, 


dé 


F is obtained: 


the following equation for 


db h e? cos(2y — 26) — e? cos(2w) —e?s*? cos(2y — 26) 
Se | ef a ee 
dt r? 15(5s? — 4) 8 
_ e’?cos(2y— 26) es’ cos(y+ 26) es’ cos(y — 26) 
60 6 2 
re 4es*’cosy 4ecosy  s*cos(20) 2es* cos(w — 4) 
3 3 2 3 
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— SS eee 
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For |2s? — 2| < 10-° 
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dt fe 15 5 15 5 
eJ sin(2w)6 a 8e cos(w — 4p) 4] 


dé h 2 - ie 
dé h 1 a, ( e cos(34 — w) - 2e cos(@ + w) * 4e cos(6 — w) ” 2 cos(26) 


15 15 as 
D. VERIFICATION OF THE RESULTS 
1. Comparison with Brenner and Latta 

The comparison of solutions obtained by Brenner and Latta must be 
restricted to the J, Je, Je*, and J? terms. Brenner and Latta limited their 
analysis in J to these terms and thus avoided the secular terms that arise in J’e’. 
No mention is made of the critical inclination in their work because |2s? — 2| does 
not appear as a divisor in the terms they included. Similarly, this present analysis 
neglected all harmonics in the potential except Jz while Brenner included some of 
the terms associated with the J, (or D, by their notation) harmonic. 

Before comparing solutions, it should be noted that there is some dif- 
ference between the two works in the notation used, this analysis preferring the 
more standard Astrodynamic notation. Brenner used the co-latitude angle (@) as 
a spherical coordinate while the latitude (6) was used here. In addition, Bren- 
ner used M as the independent variable where M + 7/2 defines the angle from 
the ascending node to the satellite. The independent variable 9, measured from 
the ascending node to the satellite, was used here. Finally, Brenner defined the 


rotation of the line of nodes in an opposite sense to that done here. In summary, 
6 — 2/2—6(where is co — latitude) 
6 —+ M+7/2(where@ is the polar angle) 
M9 -— - 


8] 


There will also be a difference in the constant terms. Brenner chose 
initial conditions that would make the analytical solutions as simple as possible. 
This analysis adopted the more general initial conditions that were outlined in 
Chapter 5. 

Now restricting the comparison to the terms that are left, Brenner’s 


solution for : and /N are respectively: 


$s ~ZJees cosy+... 


eee 


J’Me 5 4 5J 
= — 5J?Mes? + —Jec sin y + = 


2 = JMc-— 








sin(2m) +... 


A check of equations (6.35) and (6.36) shows there is an exact match of terms of 
order J. The final term in Brenner’s solution for 9 is a J* term and can be found 
in the solution for M2, equation (B.4). Next, Brenner’s abbreviated solution for 
u, 1s 


Ses? : 5 
= cos(y + 2M) + = cos(2M) ~ So ae 





These same terms are duplicated in equation (6.17). Finally, Brenner’s solution 


5 


for y; is oe — 2 which is exactly the expression obtained in this analysis. 


2. Comparison with an Independent Analysis of Polar and 
Equatorial Orbits 


Appendix A contains an independent analysis of the polar and equatorial 
orbits. Since for the polar case there is no variation of inclination, no rotation of 
the line of nodes, nor dependence of the motion on the longitude ¢, new equations 
must be derived. The equations cannot be solved in terms of the angles 1 and N. 


Instead, the equations are solved in terms of the variable A? which is related to 


dt 
dé 


For the equatorial case the inclination remains constant. The line of 


nodes is undefined since the satellite remains in the equatorial plane; however, the 


Sz 


angle M does change with time. Instead of using the polar angle @ to measure the 
angular displacement of the satellite, the angle ¢, measured from the fixed axis 4, 


is used. A new dependent coordinate is defined as 
1 ¢-K+JY,¢+ J’Yod. 


As Appendix A demonstrates, there is exact agreement between the 
special polar case and the general case for tg = 90°. In addition, when the appro- 
priate change of variables is made in the equatorial case, it agrees exactly with 


the general case at tp = 0°. 
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VII. CONCLUSIONS 
AND RECOMMENDATIONS 


The results achieved in Chapter VI represent a unique solution to the problem 
of a satellite in orbit about an oblate planet. In fact, prior to this present work 
the problem had not been solved in this representation [Ref. 18, p. 340]. In 
contrast to the universal solution presented here, all other methods require a 
reformulation of the problem in alternate variables in order to achieve solutions 
at various singularities, the most prevalent being the critical inclination. 

The results support the theory that at the critical inclination there are no 
discernable features in the perturbations of the elements that distinguish it from 
any other. This conclusion was reached by Lubowe [Ref. 38], is shared by Taff 
|Ref. 18] and has been corroborated by the physical data. However, all analyti- 
cal solutions arrived at via canonic transformations predict perturbations in the 
vicinity of the critical to be 25 times greater than those away from it. 

The perturbation method used in the analysis embodied the principles out- 
lined in Chapter I. First, the resulting solutions are significantly more accurate 
than the two-body solution. For the appropriate orbits, the relative error of two- 
body solution is 1000 times greater than the solution obtained here. Second the 
solution was obtained in parameters closely related to the classical orbital elements 
and in Cartesian coordinates; no transformation to an alternate non-physical set 
of elements was required. Therefore the physical effects are easily distinguishable 
throughout the analysis. Finally, as has been noted the solution is valid for all 


orbital parameters. 
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While the perturbation method was similar to the Method of Strained Co- 
ordinates several novel techniques for dealing with secular terms and apparent 
singularities were introduced. The verification of these techniques was possible 
due to the extensive use of the symbolic computer program MACSYMA which 


allowed for the investigation of higher order formulas. 
The analysis suggests several areas for further research. The areas include: 


e The addition of the higher harmonics of the gravitational potential, e.g., Js, 
J4, etc., to the problem. 


e The investigation into the feasibility of including non-conservative perturb- 
ing forces such as drag or solar radiation pressure. 


e Numerically integrating the equations of motion subject to the initial condi- 
tions used here in order to check the analytic results. 


e An in-depth analysis of the use of canonic transformations in satellite orbit 
prediction to include a comparison with the method used here, the purpose 
being to determine if the critical inclination singularity is an artifact of the 
transformation process. 
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APPENDIX A 
THE POLAR AND EQUATORIAL ORBITS 


A. THE POLAR ORBIT 
Reference is made to the general equations of motion, equations (5.13) - 
(5.15). Now in the case of a polar orbit, the longitude ¢ is a constant. Therefore, 


the equations of motion are reduced to 





d?r dd\* —GM Re lee 

a" (a) Soe eg (ate) ole 
d {,d0\ _GMR? | 
7 (+ 4 = —J» = sin(26) . (A.2) 


Let J = 3 (J, =>) and u = p/r 
where p=h’*/GM. 


Define a function A such that 
d6 _ 
2 =i 
dt 
so that 
d hu?Ad 
= (A.3) 
Applying (A.3) to (A.2) the independent variable can be changed from ¢ to @. 


Equation (A.2) may be rewritten as 


dA? 


“de = —2Ju sin(26). (A.4) 
Similarly equation (A.1) becomes 
d?u 1+ Ju E sin(20) + a 
Sy ee ee A.5 
do” A? 7 (45) 
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To solve (A.4) and (A.5) a solution of the form 


A* = 14+JE,+J°E, +... (A.6) 
u = l+tecosy+Ju,+J7u2+... (A.7) 
where y = 0-—w+J6y,+ J70y.4... (A.8) 


is assumed. 
Substitution of equations (A.6) and (A.7) into (A.4) and equating coefficients of 


order J results in: 


aE, 


7G sin(y + 26) + e sin(y — 26) — 2sin(26). (A.9) 


A solution of (A.9), ignoring terms of order J? or higher, is 


26 
a ecosty + 26) + e cos(y — 26) + cos(26) 


+ Kye’ cos(2y — 26). (A.10) 


Substituting (A.10), (A.6), and (A.7) into (A.5) yields to order J 


d 5e*(2y + 26 
eeu + un) = 2e€ COS Y U1 + Se"(2y + 28) = eK, cos(2y = 26) 
dé 8 
e*cos(2y — 26) _ 5ecos(y+26)  e*cos(2y) © e cos(y) 
8 3 4 
3e? cos(26) 2 cos(26) a 2 oh 
4 2 4 2 


The cos y term would produce secular terms in the solution to u,, therefore y is 


chosen as : to make the coefficient of cos y zero. 
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The resulting equation is 





2 26 2 - 
du, free Se’ cos(2y + 26) _ ?K, cos(2y — 28) + e* cos(2y — 26) 
dé 8 8 
n Se cos(y+20) — e*? cos(2y) is 3e? cos(20) is cos (26) 
3 4 4 2 
e* 1 
a Ad 
ADE (A.11) 
It is assumed that the solution to (A.11) has the following form: 
ul; = a+ a,cos(2y + 26) + a2 cos(2y — 26) 
+ ascos(y + 26) + as cos(26). (A.12) 


Equation (A.12) is substituted into (A.11). The unknown coefficients may 
then be solved for by equating the coefficients which have the same harmonics. 


The results are: 


ues ae 
0 Tie 5 
ag 
oe 
: 24 
(8K, — 1)e? 
Qo = 
8 
5e 
ee 
: 24 
e 
i = 
; 12 
e? J 
C= Bee 
4 6 


Substitution of these coefficients into (A.12) results in 


*cos(2y + 26 
= —£0s2u +2) _ 62K, cos(ay — 28) + 


e* cos(2y — 26) 
8 
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5e cos(y + 26) " e?cos(2y) e?cos(20) —cos(26) 


24 TZ 4 6 
ce: a 
Ao: 


The following complementary solutions must be added to the particular so- 
lution for u to ensure that any secular terms in the second order solution can be 


eliminated 
— Ks cos(3y — 20) — K7cos(y — 26) . 


The first order solution for u now has the following form: 


e*J cos(2y + 26) 
24 


: _ 2 6 
— ed (8K, Sb 20)  5eJ coely +26) _ J Kyeos(y — 26) 


u = 1+ecosy—J Kscos(3y — 26) — 


e*J cos(2y) 
12 


2 2 
(3e" + 2)J cos(26) et J (A.13) 
12 4 2 


+JKosiny+J Kgcos y 


The expressions for EF, and u, are substituted into (A.6) and (A.7) respec- 


tively. (A.6) and (A.7) are then substituted into (A.4). The result is 
dE e’sin(2y +40)  e?sin(2y + 26) 
do 24 12 
e* sin(2y — 26) 
12 
e*sin(2y — 46) nl 5e sin(y + 48) 
8 24 


= —K-s sin(3y Tes 46) “1 


+ eK, sin(2y — 26) + — eK, sin(2y — 48) 


in(y + 26 
— Kgsin(y + 20) + csinly + 2) 
e sin(y — 26) 


+ Kgcos(y + 20) + Kegsin(y — 26) + 5 


— Kg cos(y — 26) 


e* sin(2y) 


— Kysin(y — 46) + Kssin(3y) + e?Aysin(2y) — - 


5e sin y " e*sin(40Gz) rn sin(46) 





AP Kysin y — 4 4 6 
2 sin (20 
ie esin?h)  sin(26). 
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Integrating this expression for FE, would yield a secular term in the harmonic 


cos(2y — 26). Therefore the constant K, is chosen as —~ in order to eliminate the 


problem term. The equation is now 








dE, e?sin(2y +40)  e*sin(2y + 26) 
— —K. 3y — 48) + ———=—— — 
d6 oe 12 
5e’sin(2y—40)  5e sin(y+ 40) 
+ 7 — Kgsin(y + 28) 
26 in(y — 20 
+ esinly +2) + Kgcos(y + 26) + Kgsin(y — 20) + comly— 7) 
e sin(2y) 
— Kogcos(y — 26) — K7sin(y — 40) + Ks sin(3y) — enn 
Sesiny  e*sin(46) — sin(46) 
+ Kysiny O4 st 4 6 
* sin(26 
ar cone + sin(26). 
Integrating the above expression yields 
e*cos(2y +40)  e*cos(2y + 24) 
a 3y — 46) — eee re 
FS oe Yaa 144 a 48 
és 5e* cos(2y — 46)  e cos(y + 46) is Kg sin(y + 26) 
48 24 3 
6 
n Kgcos(y +28)  e cos(y + 26) + Kysin(y — 26) 
3 18 
— 20 K —4 K 3 
+, ceo e cos(y ) _ Kreos(y—4z) Ks cos(3y) 
2 3 3 
e” cos(2y) Se cosy e*cos(46) —_cos(46) 
ee Cee eee 
eae EOS Te 16 24 
mo cite 7 a (A.14) 


Continuing the procedure results in the following second order expression 


d’ uy eK;cos(4y— 20) e*cos(4y — 26 
eKscos(4y—46)  5etcos(4y—46)  3e*cos(3y + 46) 
4 | 576 16 
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3e* cos(3y + 26) e* cos(3y — 26) 
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12 4 
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+ See eee Ee eee 
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4 4A 32 
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It is noted that three harmonics in (A.15) would produce secular terms in 
the solution to uz. The harmonics are cos(3y — 26), cos(y — 26), and cos y. These 
harmonics may be eliminated by choosing 


3 


Ks — a 
96 
3e° — Be 
[a 
a 25e? + 34 
¥2 = ee . 


The first order solution for (A.4) is 
A?=1 + J (cos(26) + e cos(y — 26) 


ecos(y+26) | e* cos(2y — 0) (A.16) 


3 12 


u may now be rewritten as: 


e’J cos(2y+26)  5e*J cos(2y— 20) 5eJ cos(y + 28) 


= ] sO 
u +e cos y D4 a 54 34 
e°J cos(y— 20) eJ cos(y—26) e?J cos(2y)  e?J cos(26) 
24 12 12 4 
J cos(26 ees 
= = e002) — J Kysin(w — 00) + J Kg cos(w — 60) — wa = (A.17) 


J@  Je*sin(2y — 26) (25e? + 34) 
h = §- —— 4 N55 2 "7, 
where y we 5 =e 48 a 06 


Equations (A.16) and (A.17) satisfy the equations of motion to first order, 
and the equations produce no secular terms to second order. 

A comparison may now be made between the polar solution derived in this 
appendix and the general solutions for y and u, equations (6.22) and (6.37) respec- 


tively. To make the comparison easier, the relatively complex initial conditions 
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(equation 5.27) are discarded. Instead, following Brenner and Latta, the initial 
conditions are selected in order to make the expressions as simple as possible. 
In addition, constant terms like sin(2w) and cos(2w) are dropped. It will be re 
membered that these constants were added to the general solution to prevent 
singularities at the critical inclination. It is a simple exercise to add the cos(2w) 
term to the polar case, specifically to equation A.10, and then show that identical 
results are achieved with the general case at 1, = 90°. In fact, this has been done. 
However, to display these terms adds nothing to the comparison and only serves 
to make the analysis more tedious. 

Checking equations (6.35) and (6.36), it is seen for the polar case, cos 1, = 
c = 0, there results 


t=“) and Siew 1), 


This agrees with orbital theory that for the polar case there is no variation of the 
inclination nor does the line of nodes rotate. 
Setting sin 1, = s = 1 and disregarding the appropriate constants, equations 


(6.37) and (6.22) are 


e’J cos(2y + 26) Hs 5e7J cos(2y— 26)  5eJ cos(y + 26) 


a 
e°J cos(y— 20)  eJ cos(y— 28) | e’J cos(2y) _ e’J cos(26) 
24 12 12 4 
Zo 
meer 2!) _ J Kg sin(w — 69) + J Kg cos(w — 8) — aS 
6 4 2 
and 
J6 e* sin(2y — 26) 2, (25e? + 34) 
= 6§- = ———_—_— + J°6-—__—. 
y=O0-wt - ‘a ae a 96 


These results agree exactly with the (A.17) and (A.18). 
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B. THE EQUATORIAL ORBIT 
Reference is again made to the general equations of motion, (5.13) - (5.15). 
An orbit that remains precisely in the equatorial plane is possible if only even 
harmonics of the potential are considered. Now 6 = 0°, ae = 0, a = 0, and the 
equations reduce to 
d’r dd OV 


d { ,dp\ _ 


As in the general case and the polar case the independent variable is changed 


by use of the integral obtained from (A.19) 


d _ 
ad = constant =h 
dt 
d u’_ d 
— = =h—. 
7 “rr (A.20) 
Substitution of (A.20) into (A.19), results in 
d*u : 
Fa (A.21) 
2 
where: J = ra 
Z oe 
— ne 
To solve (A.21), assume a solution of the form 
u = l+ecos¥+Ju,4+J*uz+... (A222) 


where Y = ¢@-—w+J¢Y,+4+ J’¢Yo. 
Substituting (A.22) into (A.21) and equating coefficients of order J results in 


d? 
a tt =1+2J Y,e cos Y + 2Je cos Y + e’ cos’ Y. (A.23) 
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The presence of the term cos Y will cause secular terms in the solution for u, 


therefore choose Yj = —1. Equation (A.23) is now 








d*u, e? ~—e* cos? 
+t, = 1+—+t 


16? ; Sn (A.24) 


Assume a particular solution to the above equation 
ui = a, + a, cos 2Y. (A.25) 


Substitution of (A.25) into (A.24) results in the following solution for u; 


Z 


1 
uj = 1+ z - ae cos(2Y). (A.26) 


The approximation for uw is: 


a 
u=l+ecosY+J ( + - — ge COS ey) (A.27) 


The solution must be checked to ensure it causes no secular terms to second order. 


Substituting (A.26) and (A.22) into (A.21) results in 


@ 3cos(2Y 

ae te = — © S08l2¥) + et cos(2¥) + 2e cos YY, 
5e* cos Y 
WE ao 


Note that the terms with cos Y will cause secular terms in uz. Therefore, Y2 1s 


chosen as: 


The first order solution for u is equation (A.27) 
Se? 3 
where Y=¢-w-—J¢-J’¢ (5 +5). (A.28) 


Equation (A.27) is the correct solution in terms of the variable ¢. To compare 
the solution found here to the general solution for the case 6 = 0° will require that 


(A.27) and (A.28) be modified such that they are in terms of 0. 
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Referring back to Figure (5.1), it is seen that ¢ is measured from the fixed 
axis YT such that 


g=7/2+0+N. (A.29) 


Now, from equation (6.27), the first order solution for N with 1, = 0° (sin t, = 
s = 0 and cos 1, = c= 1) is 


7Je* 4 
;: sin(2y — 20) — 37¢ sin y — 








Jes ce 
+ — 


=n, Je | 
u 6° 2 


(A.30) 


(Note: As in the polar case, the constant terms cos(2w) and sin(2w) have been 

dropped in the general case for comparison purposes. As explained previously, 

eliminating these terms does not affect the validity of the comparison.) 
Substitution of (A.30) and (A.29) into (A.28) results in 


1 Jee L 


Y =@-w-—2J0+ : sin(2y — 26) — 5 





4 
J*e76 — 3° sin y. (A.3a9) 
Now, from (6.22), the solution for y with 7, = 0° is 
T 32,2 at: 
y=O0—-—w—2J0-— —J°e"6 + —Je*sin(2y — 2z). (A.32) 
12 48 
Using equation (A.32), equation (A.31) may be written as 
4 
Y =y-——Jesin y. 
3 
And therefore by use of the Taylor series expansion 
ae 
cos(Y) = cos y+ gue sin’ y. 
Equation (A.27) may now be written in terms of y(6): 
72 52 
u=1+ecosytJ (14 26? — =e cos 2y}. (A.33) 
Reference now is made to equation (6.37), the general solution for u. By 


letting sin 1, = s = O in this equation, it is easily seen that the general expression 


becomes (A.33). 
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APPENDIX B 


Znd ORDER ANALYSIS 


Substitution of (6.2) - (6.5), (6.10), (6.16), and (6.17) into (5.23) results in 


Gin dt, , aN, dt» 
2sin +79 0 8 — Gaz Sin 8 + sin te cos’? — 2-7, cos 6 


ce°K, s sin(3y — 6) 
2 
ce°K,ssin(3y — 36) 





ce°K,s* sin(3y — 6) — cKss sin(3y — 6) + 


5ce*K,s°* sin(3y — 36) 
3 6 
7ce’s* sin(2y + 36) “3 143ce’s*sin(2y+ 6)  29ce’s sin(2y + 6) 
24 72 18 
15ce*s* sin(2y — 6) 
8 
7ce*s sin(2y — 6) 
6 


+ cKss sin(3y — 36) + 


5ce’K3s° sin(2y — 8) + 10ce?K,s* sin(2y — 6) — 


4ce’K3ssin(2y — 6) — 4ce?K,ssin(2y — 6) + 


5ce’s° sin(2y — 36) 
8 
35ces* sin(y + 36) 
24 


5ce’K3s° sin(2y — 36) + + 4ce*Kgssin(2y — 36) 


4ce’K,ssin(2y — 36) — + ce Krssin(y + 6 + 2w) 


2ce°K,s° cos(2w) sin(y + 6) — ce*K,s cos(2w) sin(y + 6) 


5ces*sin(y + 6) 


2ceKys° sin(y + 6) + + 2cKgs sin(y + 8) 


4cessin(y+ 6) 


ceK2ssin(y + @) — — 2cKgscos(y + 6) 


3 
ceK7ssin(y — 6+ 2w) + 10ce Kaa coainaisioiy 7!) cos eins!) 
ce°K,scos(2w) sin(y — 6) if 10ceKos* sin(y— 8)  5ce*K,s*sin(y — 6) 
3 3 3 
co) — 2cKgssin(y — 6) — ceK7s sin(y — 6) 
ceKzssin(y—@)  ce*Kissin(y — 4) A 8cessin(y — 6) 
3 6 3 
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a (ee 
+ 2cKgscos(y — @) — ce*K,s°* sin(y — 36) + 9ces" sin(y — 38) 


4 
°K,ssin(y— 36)  5ce’s' si 3 sin(36 
+ ceK7ssin(y — 36) — ce" Ks sin(y — 38) = 5ce*s” sin(38) _ des sin(39) 
2 4 3 
5ce’s* sin 6 
— 10ce*?K,s* cos(2w) sin 6 — 10cK2s* sin 6 + EY — 
Ocs* sin 6 sin 6 
7 10cs* sin + ce” sin eee (B.1) 
3 3 
Choose Kk, = a and Ks; = or ae to eliminate the secular terms 


sin(2y — @) and sin(2y — 38) in equation (B.1), so now 





dQ d*i a di 
a2 sin im, sin 6 — ae ain 6+sin1 a cos 9 — 2—* cos 6 
39ce*s sin(3y — 0) _ce®s* sin(3y — 4) 
Sts — F555 sin(3y — 8 
1800s? — 1440 8 a 
a llce*s sin(3y — 4) ‘ 35ce°s sin(3y — 36) 5ce*®s* sin(3y — 36) 
240 1800s? — 1440 24 
7ce°s sin(3y — 30) 7ce*s* sin(2y + 36 
+ cKss sin(3y — 36) + ee omty ee ee mt yt 39) 
z 143ce*s*sin(2y+6)  29ce*s sin(2y+6)  35ces*sin(y + 36) 
72 18 24 


78ce*s cos(2w) sin(y + 6) 
1800s? — 1440 
ce*s* cos(2w)sin(y+6)  11ce*s cos(2w) sin(y + 4) 
A 120 
5ces* sin(y + 6 
* (y + 9) 


+ ceK7s sin(y+ 6+ 2w) + 


— 2ceK2s* sin(y + 6) + 2cKgs sin(y + @) 


4ces sin(y + 6) 
3 
70ce*s cos(2w) sin(y — 8) 
1800s? — 1440 
5ce*s* cos(2w) sin(y— 8)  7ce*s cos(2w) sin(y — 9) 


— ceKos sin(y + 6) — — 2cKgs cos(y + 8) 


— cekK7s sin(y — @+ 2w) — 


i IZ 72 
35ce°s sin(y— 6)  10ceK 25° sin(y — 6) ” 5ce°s* sin(y — 8) 
1800s? — 1440 3 24 
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41ces*sin(y — 6 
ateso sim 9) — 2cKgs sin(y — 0) — ceK7s sin(y — 8) 


12 
ceKys sin(y— 6) 7ce*s sin(y—6)  8ces sin(y — 6) 
3 om + ———— + 2cKgs cos(y — 8) 
39ce*s sin(y — 36) _ ce*s* sin(y—30) | 9ces°sin(y — 36) 
1800s? — 1440 8 iv 4 
See Jil 2,.3< $2: 
esa sin a88) — llce“s sin(y— 36) 5ce*s*sin(30) —_ 5cs* sin(36) 
240 4 3 
240ce’s cos(2w) sin 6 5ce*s* cos(2w) sin @ — ce?s cos(2w) sin 0 
1800s? — 1440 4 6 
Pode nuh 3: 2. = 
10cK2s° sin 6 + ee — Beeman’ _ f se . on essin@. (B.2) 


Solve for t2 and N, using the technique described in Chapter 6, Section 5. The 


second order expressions for tz and 2M, are respectively 


lice*®s cos(3y — 20)  ce®s*cos(3y— 26)  23ce*s cos(3y — 26) 


2 "900s? — 720 6 360 
_ 2cKss cos(3y — 26) — 25ce*s* cos(2y + 26) 2 29ce’s cos(2y + 26) 
3 96 144 
re lice*s cos(y— 26)  2cKossiny  22ce*s cos(2w) cos y 
900s* — 720 3 900s? — 720 
A ce’s* cos(2w) cos y _ 23ce*s cos(2w) cos y  8ceK2s° cos y 
3 180 3 
ce’s’cosy 85ces*cosy 2cKsgs cosy 2cK7s cos y 
is sc Saeco, “las 
_ 2ceK2s cosy _ 23ce*s cos y + 20ces cos y (B.3) 
3 360 9 
a, = 2ce*sin(3y — 26) _ ce®s’ sin(3y — 26) i ce* sin(3y — 26) 
75s? — 60 4 30 
4cKssin(3y — 20) 17ce*s*sin(2y+20)  29ce*sin(2y + 26) 
v 3 - 72 144 
Tces*sin(y+26) 2ce*sin(y—20)  ce*s*sin(y — 26) 
iv 36 75s? — 60 12 
_ 8ces*sin(y— 26) | 2¢eK7sin(y — 28) i llce*sin(y — 26) 
eZ 3 360 
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4ce*cos(2w) sin y ces? cos(2w) sin y ce cos(2w) sin y 


75s? — 60 2 15 


ce’s?siny 107ces’?siny 4cKgsin y 


+ 4ceK2s* sin y + ——_—— — ———— - 


6 36 3 
2cK7sin y  23ce*sin y a 28cesiny 4cKgcosy 
3 360 9 3 
5ce’*s? sin(26) in Ses’ sin(28)  5ce? cos(2w) 4 di 5ce*s? cos(2w)8 
16 12 75s? — 60 8 
ce? cos(2w) 8 5ce’s?@ 5cs?0 = ce?0_— 


be Ks76: ae eS ee B.4 
12 ie crn2e 24° 3 rie lee 


Substitution of (6.2) - (6.5), (6.10), (6.16), (6.17), (B.3), and (B.4) into equation 


(5.24) results in 


d*u» 
dg? 





120e4 cos(4y — 20)  5e*s*cos(4y — 26) 
uz = 2e cos yy, — ———_—_—————__ —- — 


36000s? — 28800 16 
7eK;s* cos(4y — 26) is 31e*s* cos(4y — 20) | 3¢K cos(4y — 26) 
2 96 
e*cos(4y — 20) e* cos(4y — 40) 9e* cos(4y — 46) 
240 15000s* — 24000s? + 9600  36000s? — 28800 
29e‘s* cos(4y — 40)  3eKss*cos(4y— 48) 923e*s* cos(4y — 48) 
192 4 5760 
eKscos(4y —40) | 89e*cos(4y— 48)  3e%s*cos(3y + 48) 
2 7200 16 
es‘ cos(3y +20) | 2le*s*cos(3y+ 20)  29e° cos(3y + 28) 
48 8 12 
3 
ote Here —soxeaey 2 
5e°s seal — 26) + 8K;.cos(3y — 26) + 17e° corse — 26) 
35e’s‘ cos(2y +40)  3eK7s? cos(2y + 26 + 2w) 
a 
3eKss* cos(2y + 20—2w)  360e4 cos(2w) cos(2y + 26) 
4 36000s? — 28800 
5e‘s‘ cos(2w) cos(2y + 26) 19e4s* cos(2w) cos(2y + 24) 


16 96 


100 





e* cos(2w) cos(2y + 26) _ 5e?K2s*cos(2y +26)  95e’s* cos(2y + 26) 


80 2 12 
5e’K2s* cos(2y +20)  65e?s*cos(2y + 26) ‘ 3eK7s* cos(2y — 26 + 2w) 
4 8 A 
3eK ss’ cos(2y —20—2w)  —_4e* cos(2w) cos(2y — 26) 
4 15000s4 — 24000s? + 9600 
172e* cos(2w) cos(2y — 26) _ e4s* cos(2w) cos(2y — 26) 
36000s? — 28800 8 
7e‘s? cos(2w) cos(2y — 26) Ms 234 cos(2w) cos(2y — 26) 
80 3600 
576e? K2 cos(2y — 26) _ 164e*cos(2y — 20) | 2400e’ cos(2y — 28) 
36000s? — 28800 36000s? — 28800 36000s? — 28800 
? Ks! cos(2y — 26) — e4s‘ cos(2y — 26) n 257e’s‘ cos(2y — 26) 
2 72 
5eK7s* cos(2y — 20) 13eKss? cos(2y— 26)  17e?K,s* cos(2y — 28) 
2 6 5 30 
29e4s? — — 20)  245¢’s? ae = Ny 2eK; cos(2y — 26) 
5eKscos(2y — 26) + e*K,cos(2y— 20)  67e4cos(2y — 26) 
3 90 3600 
e’cos(2y— 20) Ge‘ cos(2y— 40)  528e* cos(2y — 48) 

12 36000s? — 28800 36000s? — 28800 
e4s‘cos(2y— 40)  5e?stcos(2y— 46)  3eK7s* cos(2y — 48) 
——n ee ee 
GeKss* cos(2y — 40) 121e*s* cos(2y— 40) _ 23e’s* cos(2y — 48) 

4 240 240 
3eKscos(2y—40) | 29e*cos(2y— 46)  11e? cos(2y — 46) 
2 800 600 
3eK7s° col? i 2w) + eK, cos(2y + 2w) — Sp hing, conldum 2x) cosy a) 
eK; cos(2y — 2u) — 7e>s* = +46) Tes* 7 + 46) 
7K7s? cos(y + 26 + 2w) n 7Kss’ cos(y + 26 — 2w) 
3 3 
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— 2e%s*cos(y — 26) — 


240 
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es’ cos(y — 26) 
6 


2240e* cos(2w) cos(y + 20)  5e°s* cos(2w) cos(y + 26) 


36000s? — 28800 4 
7e*s* cos(2w) cos(y + 28)  7e*cos(2w) cos(y + 26) 
+ ee |, eee 
12 90 
661le*s* cos(y+20)  47es‘ cos(y + 26) 
— 10eK2s* 26) + -—————— ts 
eK 2s* cos(y + 20) + aa + ae 
10eK2s" 20)  131e*s? 20 
, eK» cos(y +26)  131e*s* cos(y+ 20) _ 5s? cos(y + 26) 
3 24 
29e° 26 
+ ae coy + 2K7s* cos(y — 20 + 2w) + 2Ks5s* cos(y — 20 — 2w) 
480e*cos(y— 26)  25e°s‘ cos(y — 26) ; 
— ee ee = 
360005? — 28800 16 cll as 


3 = 


30 
_ 1120e*cos(y— 48)  5e*s*cos(y—48) 53s‘ cos(y — 48) 
36000s? — 28800 48 12 
_ 5K7s’cos(y— 48)  Te*s?cos(y— 40)  7e*cos(y—46) —_570e* cos(4y) 
3 24 180 36000s? — 28800 

_ 5e*s* cos(4y) . 5eKs5s’ cos(4y) ‘ 19e*s* cos(4y)  5eKscos(4y) 

32 4 192 2 
n 13e*cos(4y) __1120e*cos(3y) — 269e°s*cos(3y) — 5Ks5s” cos(3y) 

80 36000s? — 28800 48 3 
i 679e*s’ cos(3y)  677e*cos(3y) | 816e*cos(2w) cos(2y) 

72 180 36000s? — 28800 
& 1le*s* cos(2w) cos(2y) _ 359e*s? cos(2w) cos(2y) 
8 240 
r 17e*cos(2w)cos(2y)  ——Ge*cos(2y)  =—=—_— 528 e* cos(2y) 
600 36000s? — 28800 36000s? — 28800 
+ Me?Kys! cos(2y) + e“s*cos(2y)  695e*s* cos(2y) _ 7eK7s* cos(2y) 
2 96 4 
_ 5eKss*cos(2y)  21e’K2s’cos(2y) _ 121e*s* cos(2y) 
4 2 240 

_ 3007e7s? cos(2y) Keene eKs cos(2y) ‘ 29e* cos(2y) 


2 800 


1137e?cos(2y)  960e%cos(2w)cos y  15e°s4 cos(2w) cos y 
200 36000s? — 28800 4 


15e°s? cos(2w) cos e* cos(2w) cos 
= — ae! aT i + 30eK2s* cos y 
4 30 
15e°stcos y 275es‘cos y . 3e%s* cos y 
16 tse ~SCOCSC AS Se ae 


6les* cos y ” Te*cos y “f 5eK7s" cos(20 + 2w) . eK5s* cos(26 + 2w) 


_—— SSS a SS a 


6 6 4 4 
eK7s* cos(20 — 2w) is 5eKe5s? cos(26 — 2w) 180e* cos (46) 
4, 4 36000s? — 28800 
5e‘s* cos(40) 65e?s*cos(46)  5s*cos(40) 5eK7s? cos(48) 
32 32 8 4 
19e*s? cos(46)  e*cos(46)  432e* cos(2w) cos(26) 
192 ~=———s«160 36000s? — 28800 
1056e? cos(2w) cos(20) 3e*s* cos(2w) cos(26) 
36000s? — 28800 8 
e*s* cos(2w) cos(20) = 19e*s? cos(2w) cos(26) 
2 80 
23e7s? cos(2w)cos(20)  3e*cos(2w) cos(26)  11e? cos(2w) cos(26) 
120 200 300 


276e* cos(26) 7e*s* cos(26) 
emit 19d Bono nls anes Cos(00) 
ems? sansnn a 2 ce) pence 28) 16 


127e*s‘ cos(26)  7s*cos(26) 9eK7s?cos(20)  3e?K2s’ cos(26) 
16 6 2 2 
193e4s? cos(26) 95e7s? cos(20) 2s? cos(26) 
Rec OO AAR AS el Ra Sa 
+ K2s° cos(26) + AG 5 - 
19e* cos(26) 2e* cos(2w)? _ _50e* cos(2w)? 
300 15000s* — 24000s?+ 9600  36000s? — 28800 

e*s* cos(2w)? 23e4s? cos(2w)? Te*cos(2w)? 576e?K, cos(2w) 

32 960 3600 36000s? — 28800 


— 4eK7cos(26@) + 


304e* cos(2w) 2400e* cos(2w)  e*? Kys* cos(2w) 1 17e4s* cos(2w) 


36000s? — 28800 36000s? — 28800 Zz 24 
Ses‘ cos(2w)  3eK7s?cos(2w)  3eK 5s? cos(2w) 


SSE OOOO i OEE === — 


Z iz 2 
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19¢e?K»s? cos(2w 533e1s? cos(2w 41e*s? cos(2 











60 720 24 
+ eKscos(2w) — eo) st se cost?) _¢ cost) 
T SO ™ at + 2K3s* + — + 20K. s' 
+ — a ost 4 28, ME Be a, (B.5) 
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APPENDIX C 


THE CRITICAL INCLINATION 


This appendix contains a discussion of the critical inclination to include its 
mathematical and physical basis and the methods used by various authors to deal 
with the problem. 

In Chapter 6 it was found that the divisor ( Ssint a= ) often appeared in sec- 
ond order terms. A brief explanation was given for the appearance of this divisor, 
and it was shown that the apparent singularity may be dealt with in this analy- 
sis in a straight-forward fashion. However, in other methods the problem divisor 
has caused enormous complications. This divisor or “critical inclination” problem 
appears to be universal. Hughes’ investigations (Ref. 32] reveal that every ana- 
lytical theory on orbit perturbation analysis contains the critical inclination as an 
inherent characteristic. As stated in Chapter 2, the divisor causes an irremovable 
singularity in Kozai’s method. Jts presence has prompted many authors to devise 
very elaborate techniques to get around the problem with the result that many 
otherwise elegant methods are somewhat diminished. 

The mathematical reason for the appearance of the problem term is as fol- 
lows. It is recalled that one of the secular effects of the Earth’s equatorial bulge is 
to cause the line of apsides to rotate. A first order approximation for this rotation 
is given by the rate of change of the argument of perigee which from Allan [Ref. 
38) is 

dw 3Jon R? 


Bio Pics 
ae aT eee (5 cos: — 1) 


It is readily seen that if : = 63°26’, the perigee position is stationary and the 


line of apsides does not rotate. When the usual Poisson method of successive 


approximations in terms of a small parameter is adopted for perturbation analysis, 
the above equation appears in the denominator of higher order terms. 

Coffey [Ref. 34] gives an excellent background discussion on the history of 
the critical inclination problem. He states that A. A. Orlov (1957) was perhaps 
the first to notice the unusual situation at the critical inclination. Orlov found 
that the conventional Lindstedt-Poincaré algorithm failed at this critical inclina- 
tion and at zero inclination. Nevertheless, Orlov ignored these exceptional cases 
and developed his theories without addressing the problem of singularities. Krause 
(1952), Roberson (1957), and Herget and Musen (1958) had all overlooked this 
difficulty in their assessment of the long term effects of the Jz term on the Kep- 
lerian ellipse before Brouwer (1958) finally pointed out the problem to the latter 
authors. Brouwer hoped that his use of von Ziepel’s method of eliminating vari- 
ables through canonic transformations would allow him to dispose of the terms 
that lead to singularities. However, as was pointed out in Chapter 2, his solution 
contained singularities at, not only the critical inclination, but also zero inclination 
or eccentricity. 

Finding no method to fit the critical inclination problem into current theories, 
astronomers have devised alternative theories to deal with orbital inclinations near 
63°.4. The critical inclination, they reasoned, was a small divisor problem and 
therefore could be handled in the standard way introduced by Bohlen [Ref. 35]. 
That is, in the vicinity of the critical inclination, the solution may be expressed 
in terms of the square root of a small parameter, i.e., VJ instead of J. [Ref. 36]. 
The resulting solutions are totally different in character to the non-singular case 
[Ref. 30]. For example, instead of precessing secularly the argument of perigee 
will tend to oscillate about a central position. Taff notes the lack of mathematical 


rigor in this method of handling the critical inclination by stating that expansions 
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in VJ and theories about the libration of perigee are a result of misapplication of 
perturbation theories (Ref. 18, p 340]. 

Usually when singularities arise in celestial mechanics there is a physical 
reason for their occurrence. For example, singularities in the solutions for the 
motion of asteroids and natural satellites are caused by the fact that their orbital 
periods are nearly commensurate with that of the disturbing body. Under the same 
circumstances, an artificial satellite may encounter resonant perturbations in its 
orbital elements caused by tesseral harmonics. However, such a physical reason 
does not seem to suggest itself for the critical inclination problem. Message (Ref. 
37| implied that there may be a resonance between the satellite’s mean motion 
relative to perigee and its mean motion relative to the ascending node; however, 
most skeptics have not been convinced [Ref. 32]. If there were resonance one would 
expect to be able to experimentally measure deviations caused by resonances, 
but satellite guidance engineers have reported no adverse effects at the critical 
inclination. 

Another check on the physical effects of resonance should be a numerical 
check of the equations of motion at the critical inclination; however, this too 
has been inconclusive. Lubowe [Ref. 38] carried out a study in which he inte- 
grated a satellite’s equations of motion for periods up to 24 hours with various 
initial conditions at inclinations near and away from the critical. He found no 
discernible features in the perturbations which distinguished the critical from any 
other. However, in a separate analysis, Hughes came to the opposite conclusion. 
Hughes developed the Hamiltonian in terms of the Hill variables and showed that 
the critical inclination remained an inherent part of the solution throughout a 
multitude of various canonic transformations. He then integrated the equations 


of motion numerically and was able to show some resonant effects in the perigee 
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height for satellites at the critical inclination. He then reasoned from this analysis 
that the critical inclination is unique, that it represents actual physical resonance 
and that it is not merely a by-product of the method of solution or the type of 
variable used in the analysis. 

This disagreement characterizes the lingering debate over the question of 
whether the critical inclination is an artifact of the analysis process or a physical 
reality. Most textbooks choose to mention the problem; however, few voice an 
opinion. Roy [Ref. 20] devotes only two sentences to the subject, merely noting 
that some perturbation analyses break down there. Hagihara [Ref. 39], mentions 
that the critical inclination is an essential feature of the perturbation process 
but notes that heated discussions continue concerning its physical reality. Other 
authors, namely Herrick (Ref. 40], Kovalevsky [Ref. 41], and Geyling [Ref. 1] 
mention the problem and the standard procedure for dealing with it, but they 
make no mention of physical effects. Taff [Ref. 18, p 340] takes a firmer position 
by stating that it is not a prediction of second-order perturbation theory that 
there are infinitely large or infinitely rapid changes in the orbital elements, and 
that the correct resolution of these “unphysical” predictions is not to rely on bad 
mathematics. 

This present analysis lends support to Taff’s assertion. Although the criti- 
cal inclination problem was manifested in the perturbation process, it was shown 
that it does not cause singularities, nor are there any unusual physical effects pre- 
dicted by this theory. In fact, it was a remarkable aspect of this analysis that 
in the limit as the inclination approached the critical, all potential singularities 
were canceled. This is not true of the analyses which use the technique of De- 
launey/von Ziepel canonic transformations. A remarkable aspect of the canonic 


transformation method is that the non-integrability of the equations of motion 
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is not obvious when the system is expressed in coordinates or elements, but after 
the canonical transformation it is quite clear [Ref. 42]. The critical inclination is 


an intrinsic singularity of the method [Ref. 32]. 


109 


10. 


ile 


i 


LIST OF REFERENCES 


. Geyling, F. T. and Westerman, R. H., Introduction to Orbital Mechanics, 


Addison-Wesley, 1971. 


. National Aeronautics and Space Administration Report R-99, Theory of Sec- 


ular Vartations tn the Orbit of a Satellite of an Oblate Planet, by W. A. 
Mersman, p. 34, 1961. 


. Brouwer, D., “Solution of the Problem of Artificial Satellite Theory Without 


Drag,” The Astronomical Journal, v. 64, pp. 378-397, 1959. 


. Aeronautics Report U-2556, General Perturbation Techniques for Satellite 


Orbit Predtctton, by J. L. Arsenault, J. D. Enright, and C. Purcell, v. 2, p. 
6, 1964. 


. Smith, O. K., “Computations of Coordinates from Brouwer’s Solution of 


the Artificial Satellite Problem,” Astronomical Journal, v. 66, pp. 359-360, 
1961. 


. Lyddane, R. H., “Small Eccentricities or Inclinations in the Brouwer Theory 


of the Artificial Satellite,” Astronomical Journal, v. 66, pp. 555-558, 1963. 


. Sterne, T. E., “The Gravitational Orbit of a Satellite of an Oblate Planet,” 


Astronomical Journal, v. 63, pp. 28-40, 1959. 


. Garfinkel, B., “The Orbit of a Satellite of an Oblate Planet,” Astronomical 


Journal, v. 64, pp. 353-367, 1958. 


. Garfinkel, B., “On the Motion of a Satellite in the Vicinity of the Critical 


Inclination,” Astronomical Journal, v. 65, pp. 624-627, 1960. 


Vinti, J. P., “New Method of Solution for Unretarded Satellite Orbits,” 
Journal of Research of the National Bureau of Standards B, v. 63, pp. 169- 
204, 1961. 


Vinti, J. P., “Theory of an Accurate Intermediary Orbit for Satellite Astron- 
omy,” Journal of Research of the Nattonal Bureau of Standards B,v. 65, pp. 
169-204, 1961. 


Kozai, Y., “Second Order Solution of Artificial Satellite Theory Without Air 
Drag,” Astronomical Journal, v. 67, pp. 446-461, 1962. 


110 


7 eS =_Z«  —ee.gXKbea_ Oe 


13. 


14. 


15. 


16. 


17. 


PS: 


19. 


20. 


21. 


22. 


23. 


24. 


20. 


26. 


27. 


Morrision, J. A., “Comparison of the Method of Averaging and the Two 
Variable Expansion Procedure,” SIAM Review, v. 8, pp. 66-85, Jan. 1966. 


Liu, J. J. F., “Satellite Motion About an Oblate Earth,” AJAA Journal, v. 
12, pp. 1511-1516, 1974. 


Jet Propulsion Laboratory Report 32-1513, Methods of Averages Expansions 
for Artsficsal Satellite Applications, by J. Lorell and A. Liu, 1971. 


Computer Sciences Corporation Technical Report TR-77/6010, by W. D. 
McClain, 1977. 


Hoots, F. R., “Theory of the Motion of an Artificial Earth Satellite,” Celes- 
tral Mechanics, v. 23, pp. 307-363, 1981. 


Taff, L. G., Celestial Mechanics, John Wiley and Sons, 1985. 


Arnold, V.I., Mathematical Methods of Classical Mechanics, Springer-Verlag, 
1978. 


Roy, A. E., Orbital Motion, pp. 292-294, Adam Hilger Ltd., 1982. 
Baker, R. M. L., Astrodynamtcs, Academic Press, 1967. 


Kozai, Y., “The Motion of a Close Earth Satellite,” Astronomical Journal, 
v. 64, pp. 367-377, 1959. 


Handbook of Astronautical Engineering, pp. 8.2-8.37, McGraw-Hill, 1961. 


Danby, J. M. A., Fundamentals of Celestial Mechanics, The Macmillan Com- 
pany, 1962. 


Roberson, R. E., “Orbital Behavior of Earth Satellites,” Franklin Institute, 
v. 264, pp. 181-201, 1957. 


King-Hele, D. G., “The Effect of the Earth’s Oblateness on the Orbit of a 
Near Satellite,” Proceedings of the Royal Society, A, v. 247, pp. 49-72, 1958. 


Brenner, J. L., and Latta, G. E., “The Theory of Satellite Orbits Based on 
a New Coordinate System,” Proceedings of the Royal Society, A, v. 258, pp. 
470-485, 1960. 


28. 


29. 


30. 


31. 


32. 


33. 


34. 


35. 


36. 


37. 


38. 


39. 


40. 


41: 


Nelson, W. C. and Loft, E. E., Space Mechanics, Prentice-Hall, 1962. 


National Aeronautics and Space Administration Report TT F-579, Motion 
of an Arttficral Satellite in an Eccentric Gravitation Field, by V. G. Demin, 
p. 89, 1970. 


Milani, A., Nobili, A., Farinella, P., Non-Gravitational Perturbations and 
Satellite Geodesy, Adam Hilger, 1987. 


Stanford Research Institute Report 1s An Approximation of Polar 
Orbits of Near Satellites Around an Oblate Planet, by M. Weisfield, 1960. 


Hughes, S., “The ‘Critical Inclination’: Another Look,” Celestsal Mechanics, 
v. 25, pp. 235-265, 1981. 


Allan, R. R., “The Critical Inclination: A Simple Treatment,” Celestzal 
Mechanics, v. 2, pp. 121-122, 1970. 


Coffey, S. L., Deprit, A., and Miller, B. R., “The Critical Inclination in 
Artificial Satellite Theory,” Celestral Mechanics, v. 23, pp. 365-405, 1985. 


Jupp, Alan H., “The Critical Inclination Problem with Small Eccentricity,” 
Celesttal Mechanics, v. 11, pp. 361- 393, 1980. 


Hori, G., “The Motion of an Artificial Satellite in the Vicinity of the Critical 
Inclination,” Astronomical Journal, v. 65, pp. 291-300, 1960. 


Message, P. J., Hori, G., Garfinkel, B., “On the Critical Case in the Motion 
of an Artificial Satellite,” The Observatory, pp. 168-170, 1962. 


Lubowe, A. G., “How Critical is the Critical Inclination?,” Celestial Me- 
chanics, v. 1, pp. 6-10, 1969. 


Hagihara, Y., Celestral Mechanics, p. 431, The MIT Press, 1972. 


Herrick, S., Astrodynamics, v. 2, pp. 277-279, Van Nostrand Reinhold 
Company, 1972. 


Kovalevsky, J., Introduction to Celestral Mechanics, pp. 91-94, Springer- 
Verlag, 1963. 


PZ 


42. Ferrer, S. and Williams, “Simplification Toward Integrability of Perturbed 
Keplerian Systems,” Annals of the New York Academy of Science, v. 536, 
pp. 127-139, 1988. 





INITIAL DISTRIBUTION LIST 


No. Copies 


Defense Technical Information Center 2 
Cameron Station 
Alexandria, Virginia 22304-6145 


Library, Code 0142 2 
Naval Postgraduate School 
Monterey, California 93943-5002 


Commander 1 
Naval Space Command 

Attn: Code N155 

Dahlgren, Virginia 22448 


United States Space Command 1 
Attn: Technical Library 
Peterson AFB, Colorado 80914 


Director A 
Navy Space Systems Division (OP-943) 
Washington, DC 20350-2000 


Superintendent 1 
Space Systems Academic Group, Code 72 

Naval Postgraduate School 

Monterey, California 93943-5000 


Dr. Claus Osterwinter 1 
Naval Surface Warfare Center/K104 
Dahlgren, Virginia 22448 


Professor Donald A. Danielson 2 
Naval Postgraduate School, Code 53 
Monterey, California 93943-5000 


114 


10. 


11. 


Lieutenant Christopher Sagovac 
Naval Postgraduate School, Code 53 
Monterey, California 93943-5000 


LTC James R. Snider 
P.O. Box 961 
DeQuincy, Louisiana 70633 


Department of Aeronautics and Astronautics 


Naval Postgraduate School, Code 67 
Monterey, California 93943-5000 


115 


















a LS led Oe iti beail 
emp TTT Leh teeny Pe. sdlhrdn hd lant ag 
rr rh ald PPT lad alee ry Aah ; - 


Pee PP eT Ed eee yeh ed oF. a 
ee ro ad Peerar ea Pt ee eet ad 


PY lL a api BAcees 




















Satellite motion around an oblate planet 






















et eat OE OE DORE rh PL 
OTS O TY PO taka Te eater 
yh ok eel ny td eadide! EL cee PERTH 
OI ae tee al De ater 
ae ene PITY i hited bats DOLRUOU Le 
 ehaalh af nd aeiinamie Peer Baden PO A OT lal tel et ed ' 
PTT TTT Lan became ba Pere ree ae PIPER ho eee ; i 

pee PT dL Ad rT aofes over scomontter ss TT ae , 

oor ret Ak PPT Le aerrreeer ny hie Thies a pat é / 
ee ee parr wr eT Pld ded Ory TT Lael ae rere tl Tet chen : 
et ld er Te es Mere? matgtet 4 Gon vet ; 

PP Reel hie cra PPT A Toke yearn? j 
eee cee? e Py ats le rd | 
| 










pt ae Bye 3 2768 000 83662 1 


Prey ed 
yee 


Pra 
* 
Pee Mead 


aye Part Cir ce ee : DUDLEY KNOX LIBRARY 


- i 
rr 


O 
5 =Os 
i. td 
’ i) 









pet dd 
ye ad er ed oa 
TTT deen hl ke dated 
eT hd Fo moe 

et er kel bn epcobels 
oe a ad Le Foe 
gag a Borer pages? nae aT Dhak 





Lee el 





aoe 
Powys Fafear- ogi are® “ Nd 


PITT 1 iad di CPT et id ed el 
baer PPT CE 2) LR and aes Perro 5 Hoesrete® 
’ 


POET TO eh aihrte nO let 
Pr Perr 
























































































































































































































































ee a 
"sella eee rT Li cad etd te 
TT ta dnad rd PT Tt od ah eh 
dl OP at Ae Le on aw 
~ Pal ef os Peres pacet Octo! 2% oF F F Pi ec 
PPP rT ti ee ad Py fae Aaa 2 : 
PP el, eo atetenntrrze ° eet r n ’ y 
ae AT ddd Pee ee o on emcees 3 re P 4 
papacy i TT kaa pepeggceet Ces rr ae Ae ae i ; : F : ‘ , 
pale Per eat be | sane, SOONG eae PP : , 
epee er sree priya Cid data ae 6. Arn delle te Dh HR ; ; “ 
Pe Et ke Serre ih wnzee eet? re ae A Pree sane - sare | A A cE A a) ' 
ied Le PPP ade ee er er Ti PPC of pt coco scenes? eg * aT a A 3 nee a eis er P A A A ¢ 
Te er es Td Parry rire Tr ee sl Went hd ae hr te Mk ope P Atay ne Ta os 5 y 
at dteal ert ald Peodoop bor eme 0-89 %8 ey a Pe rr . : A ' a 
ere nad PRP Pt ad pO dl wm lie PTT ts das re 5 Tete ® 
err tt ETT, Tata art Pee ad fe) had Th 1] ee ee iu 
Te Lateaal yt ed A deel shore gad Pd od hel Lar ote . ao 
a Ary ee ee vos © ged a ee Wa ’ 
OP OPP PO I aN ter PT ee pee Ca o e ' 
apa) ane rr eee Wt a Pon at F 5 y 
er rd oted ad er Ya hein , s 5 
‘me 0 Cer e* TL ah A ery ra Py ere) Pad ‘ r i ) . a 
oF gh Gt we a0 ae Dah eh id OF ateTr Pr dead Arye is all As FP , Py ' 
a oma a ron 27 0h Pe eore 5 ee) Pru O a) I 
Pt eta OF et hed Pa e 5 wee iY e ‘ Oo 
ee Te te etd Wats Peer eaes ve ® E Re a ee ry ‘ og * 
ed ot dete of 7 PTT cde o A f te ¢ UJ 
ign I oT eral ’ a ed ee ' na ‘ 
Pe ed be a sets ears A re ee .? 8 ' ‘4 es 
pa eed ee Td week aide Pt ory ee og egmreet o:t 1 fn rear | ‘i 
PT Na td od di ded 5 er re alk Pree a aad ert 16 n D A 5 
ae al PO eed 7 ATE Al ie PP de a ’ ' i LU ‘ ‘ 
PIO Pr ad eid OP A f e aT) Cd ‘ Pi - ‘ ° ° a 
owe ae P Pd oe) Por fu ponders” ° ria) e Pied ene Fi f 
Te had del Py hed ed PU ed Lotta alah ae - PY abel vs 5 n soe ’ o o 
ant ed eal PP PTT Roied fete oO nase! 5 , 
A A iS ad 
i ee aioe Pr tee Sere POEL dar" Dibdnee aa F 
eed ee Py pee 1 erry Se c hi rn ' ‘ u 
«meee? was gtesey of # tem treet Pe Td ; . ° t of ‘ | D ' 
— ee ee ae] seybon op seyeqnggr? Cd a r ' . 1 1 
get re To ele) eee Pe rer Pd ae ; te ates ae . 8 f ' 
are Py at dled pele, ee aU ha Pr fit) fe 7) aad eo Us p f =a O ' 
a Tal tailed A a oe ee Pte ae al er Pe ed Fy er ‘ ‘i o Py 
Pe tad eer eT om: 1008 . oreete a aesaeqar” * eo re 7 rar ad ' r 
I Acad ee ae nde Pe iy a = aa obs erat ds ‘ oa = r n ‘ n : 5 
eee dl Jt dt ded AY Cry Lt tl il Peers ‘ee y = i A a , Par - 7 
aa Pe Oe ed A dh © yargee ote? tae we weoter * id Me . ’ : 1 ee? A : 
4 * ie ee ed I a atid ‘ ae e. 3 ‘ 1 A se 1 , 
Pe he el er eee ledae 097 Fe? erase : Ld U Sr ; i ss 
Fs - pl ipe te PA ee Rd Fs ene Pn ae) 
°° Pd ROE FOr er i cd PT 
= > és = Fy ase ‘ rr ef oa ' ' oe 
er or e Pr re ee oe a U . 
a ad DY ed id adale Perry 1 beale A e 7 r ee eae 
cad te tar Ff we rea Pee ae dad ere gs ate ft rv ; n ‘ ml , ; oa] f 
Sure ome Pe a dt PP a Us te ee uote wh 7 , oe Zs A F 
- A A P , 
al oe ato 78 FoF Pe Pa Pd ' 7 ‘ 4 ar 
* TO A Sool odie - ‘ an eo? oe . ‘ 
. . a iu Pan Pat 
a on 8 ter @ ie 
=. #F« - be s : , 7 
5 he a ? n by f . y 
ape vs A A ” t S a 
ea Plow niet io 8 5 U el ee . 
Pa d + A A . $ U 
n f A a Pa | of ° F 
‘ Py Cd ' LU é es 6 1 sd 
; , r , . ‘ : 
ee p a I 5 Pi ' 
a , 
o oan dg Gofal 8 - ° LE od i 
y eT ee a Per) p P r 5 ‘ ® 1 Ps 
Fy ot ot ae A ed Pare ’ or 2 ‘ a ‘ . r} ‘ 
e rd Pry SS Pi i a. ie : . is 5 ° : LY : 
- b On o# ig Pru mere iy 1900 o e- a es ao : ‘ Z : 7 7 
el 7 . 2% Fy : ‘ A re Ar ie Set) ° ' tes O A 
des pa Ch aa A Arr Pd he soy : f - " , H 
ear o o . Py) ° - A a e? 5 : ‘ ' ° ran) 
- = a ee HS OE emit © ee ¢ Uy n 5 
= ¢ * ‘ 7 rte! 0 3 ar Pn . cf Lan) 
“ag > i ee tial * nen A ra oa es ‘i r A 
Le ” ir ar o- ° Pay : ig f ee 4 ae =e 
Oe? pe ore f ted are : : ‘ K : ie : : 
el - ry oa Pa. ye Peer ‘ Ls ates ion : g : J 
aT ® Pe 7 . ° ms ‘ 6 r .” - + ehpet : . is ss A 
- = oo F hen ee ree a rn Ce ; “ 
= eee A i n : . ; 
d et ore cd a5 oe , © . CULL iD ‘ fee 
i : sf a) > ; oe ” ? ; ? ° e n 5 ey ae 
a4 pet F] H s° ° ' ’ Pa - Pen . n o of a 
Aw. F = yes A U ' : Le ‘ an : : 
oO ‘ a Be Ld | A A A ry oT : ' . 
ra Fi " | mers ° Fret e 2 rr Pry) Py Py rim 1 or fe ‘ 
iG i. A Pa : | ‘i 5 é ‘ D 
- , tL et tte Ta Ad =¢ a U a s 1 
o Pe Se | ee aa i els Ui . 5 ne o Ed ee i 
7 # ° A core F : a] eh ewe ad o* “ 
att a ee ay re er eae = ry o * , D a 
-" o = -_) re a bet ” ’ A o* Ps a F rr . A A A PS 
; a ° Pap caer i had iP et ° P 1? 5 A A 
A ara Lae re - » y U ait ek Jog Oy ey aha 
= ° PF ea ae ° Hi P . F f APH Lies hd a C ‘ Au ° : 
2 s a ae ape ¢ : a ar oe . 4 f ; ce a ; he 
ba page” ee aot ps ° . r . : A ‘ . o A a 7 ’ 
ad a ace Fa ry a = rt > Je ° 7 fi . D a 4 7 
ae Fy se P AA r zi A Feri. ree é - = ° fe . . 7 
« oT) ¥ ae = . re et ae F Pir Pa ee A 1 ° . 1 
nic , p Pr ” 3 oe - SC te tae i eee Ui : ‘ 
+ » or 3 2 ee he " . y Lg eee: ‘ . 
= F a Pe al ted Perce iar on) oe a re ee Py A o 4 ee oe a 
Ces ? z Ul “r°# 4 eT) %. ¢ = = ‘ i Fi 
-* a ‘, ry = j a Ps ef ry Pa nd ° er 8 Fi 
Bs ore wef i La Pe eee) P vt or a F Oo : ° 
° r ’ P S F 
2 oe EI Ta Ps . ° Pad ° aa ‘ , si 
° = a F j Yeo Per : S ° nome | s , : 
- . ° e a r 2 bt" at Fa es : 
3 ? oy of “ Pde os Aa . A 5 D F 7 
Cs P 5 oad a Ps rf - PP us ms 
Cd a eo eee 8 a ye od Fj Q c 2 n - é- . Py (wa U A 
i 2 j : Pin ae ed Boy uC 7 4 . 5 , 
2 4 = (oe Py P aan 1 P Pa oe A an 
* . F $ P . P D oO . 
a7 Peary i ° o * FP Pe re Lie, 5 p ‘ ' r P é s 
ae ae - PY e * eta) P ° = = 7 O 
= , ps ai ae 5 P a oa “3 ' ~ . 7 
oo * , 2 r , S : Ss 
Ca Py - - @ . . . zs J Ul 4 hr ee ee n D a 
a = ry Bs asa e 8 ae PF y 7? 1 ¢ ' 
€ a r , = F r A td ° % ' ° p ) rf : ' i 
Pi 3 = S ’ Pas LU 
bs a ° fT A Ps 0 ‘ r 
: a 2 P er » co r Ps . r Cf 
F # >. » * 6 A en te Pi L . 5 . o Fl F 
Py . r ee 5 
» P s ) ar Fi A FT 
a oe 5 - 
® : a , 5 
= 7 . D 7 r ea f o 
, a 
> > wae © iw - rae 7 
7 ss od H . U ad J . " 
5 3 F we 
- ’ 
a he " . . - 
a * od + ; 
Pe 
y , oF ji 
- q a _? P ‘ , 5 
7 eo = . ? p 
] ° o a at A es § 
« rh ~ ° 4 
° Co ca a rf . J ® id 
, P t 
, i A F i he 
a = bel dd s ‘ Co « 
P sf J ° f ' 
Cn = 5 = PP ‘ Fa A Ey 
rs *-, Fa 23 F Pa P @ ‘ ' . ® 
° Mey cs 7 p O 
" : f * : Tt ‘ t ‘ ° 
* & . P s oi i P ae ' o 
° 
* ’ bf + . H] 
° oe A = 4 - 
a5 P - 5 ra . 
= al a td A ] ad , 7 ® | FI D ‘ o 
= i «4 * > u . in 
C ie Tr ss = Fy 
* ° = s 4 *s wr fi ‘ ? Fi id * bs 
. a « A * a ry ° . : 
dled , P S , A pi , . ‘ 
~ . ° - F : on e¢ o ¥ ° n 
: 6 B , reer F 
° a ~~ r ® bd ri ° ¢ : 
od Pa i Pe ee r ' I . 
a ee se ° rs e Y P . A 
Lf F *. e c ° 
cd i. ba , Ad bd ® 
« ‘ cd Pa at 2 ly : 
: al « on PT The 7 A "y . E « , ° 
¥ . “ w @? a 3. s ad a ‘ c ’ = Z : 
ad ° ™ n a A TT € bo ® rc ' 
| ° J . Y 
q . ’ ~ - C7) ’ 4 7 ‘ 
¥ e ‘ - Cr «= L € ‘ ® ¢ 
. . - Py ? S . 
Ps - er ee ' 
or , my io ofe - Ls 
- i S . ‘ | % < ‘ Y s 
- a Vote Bid o 7 La ry ’ t 
, o 7 r ° Pa | s Y bd bd 
Ps e rs “ "i « ae? 4 ‘ ' Py 
P es . ri F a] mL ha 5 
= - ve e ° D . F ‘ 
le , a e iP o fs 7 p 
. . " a) F f Y J 
e ® . ° » 8 i 1% he] ’ 2 $ 
r , °  ) f ;? P f p 
» y 5 oP . me, oe A Y c 
ms Py e " P iu ‘ 
* J La ad 4 fs t 
‘L . : 5 
v P P 
° 
J ¥ s & Ff . f , : 
YT) es 7 U a Py ? ) ‘ 
a - e 
A f F a 
rh) ad « 
. hd . r ® ’ vr eo:@ ‘ A 
» ° ~ t CU ry r ri ry ir Ln) " ‘ ¢ . ‘ i fa es be Mi 
a oe] , Py" e Et Pie) 4 kL 5 y 
e r es > a  # , ci ° i] 
ie Ad . e °alf . ‘ Cat 
a © . CY s A 5 ¢ e 
r a eT r » ¢ 4 Ae P pI 
P = 
- . D t « S 4 
bd ‘ be . 7 « ® Y be ad a ’ ‘ i ‘ 
2 Pr ¥ a f opp iwreem A i ‘ ‘ 
e si t wv ° e* j -*® LI L cf 
O ‘ e n F ’ 5 
bs 6 "hed * f ? P : ‘ ? 
: od] a » r P eoctes 8 r 
pair! ¥ 2 Fe . ca e’Wv a e rt ] "| is 
aol tees | oy ry r a) ry 7 
My a healed ve -7@ Ps r " . e he 
Sty a ) c ao « 2 > ’ . ' ‘ e 
y hy a] Ta Ce a 4 " Pe tors ’ f 
. =e Stig . ® ro ry P oy . a ‘ U 
; » To © ] oe f 7 Y ? r ty cs 
- m ed 4] A - rn , aa J Py @e . 
5 se - “a t ‘ rrr iL aa ran) L 
eer ® . Py « “ Cy e C e ' 
-, ted . e 4 . vwe!l eet ’ ed r 0 ' e 
. Pan a io Fi ‘ ry fn ‘ t 
e . A f a ott an) « ’ 
r 7 ee ri wr) ‘ * 
L tacheabintin c S * e oy cy A 
= i ¢¢ re ro | P i i ‘ 
a PT r ‘ 
nae be af ? d v 5 
ie pe Owe P) - ee * 
4 e , 4 is ad ¥ oe hae 5 
P ns dhe pp bed . P " , i] 4 Fi 
ae ote 4 P f ‘ ‘ 
C viii Pre? ot) PY o 
ag ] ri i) : - 
: . 5 
Ld ~ e a a 
J 2 . ‘ P , 
; é ‘ , . ~ 
into tied e bs bday bad a se é e* f LY 
Bode Ad rr % oT a as 
ed r rar) , + re 
CG 
Py 5 
we rs " PF ig 
- , Py . - ps P ’ 
hai Del . ad a Ch ] ov | P rE ‘ 
r rr ee 5 f 
r aad a ’ 
Dd ¥ Tee Ad eae ° if n 
a e ¢ a 2 
i be . 
‘ 


